source: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc @ 1312

Last change on this file since 1312 was 1228, checked in by garnier, 15 years ago

update geant4.9.3 tag
  • Property svn:executable set to *
File size: 25.7 KB
Line 
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,v 1.20 2009/12/16 17:50:11 gunter Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33// File name:     G4GoudsmitSaundersonMscModel
34//
35// Author:        Omrane Kadri
36//
37// Creation date: 20.02.2009
38//
39// Modifications:
40// 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
41//
42// 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta
43//                     sampling from SampleCosineTheta() which means the splitting
44//                     step into two sub-steps occur only for msc regime
45//
46// 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
47//                     adding a theta min limit due to screening effect of the atomic nucleus
48// 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
49//                     within CalculateIntegrals method
50// 05.10.2009 O.Kadri: tuning small angle theta distributions
51//                     assuming the case of lambdan<1 as single scattering regime
52//                     tuning theta sampling for theta below the screening angle
53//
54//
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56//REFERENCES:
57//Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
58//Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
59//Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
60//Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
61//Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
62//Ref.6:G4UrbanMscModel G4_v9.1Ref09;
63//Ref.7:G4eCoulombScatteringModel G4_v9.1Ref09.
64//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
65
66#include "G4GoudsmitSaundersonMscModel.hh"
67#include "G4GoudsmitSaundersonTable.hh"
68
69#include "G4ParticleChangeForMSC.hh"
70#include "G4MaterialCutsCouple.hh"
71#include "G4DynamicParticle.hh"
72#include "G4Electron.hh"
73#include "G4Positron.hh"
74
75#include "G4LossTableManager.hh"
76#include "G4Track.hh"
77#include "G4PhysicsTable.hh"
78#include "Randomize.hh"
79#include "G4Poisson.hh"
80
81using namespace std;
82
83G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
84G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
85G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
86G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
87G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
91  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV),isInitialized(false)
92{ 
93  fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
94  particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
95  tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
96  tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
97  theManager=G4LossTableManager::Instance();
98  inside=false;insideskin=false;
99  samplez=false;
100
101  GSTable = new G4GoudsmitSaundersonTable();
102
103  if(ener[0] < 0.0){ 
104    G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
105    LoadELSEPAXSections();
106  }
107}
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
109G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
110{
111  delete GSTable;
112}
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
115                                              const G4DataVector&)
116{ 
117  skindepth=skin*stepmin;
118  SetParticle(p);
119  if(isInitialized) return;
120  fParticleChange = GetParticleChangeForMSC();
121  InitialiseSafetyHelper();
122  isInitialized=true;
123}
124//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
125
126G4double
127G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
128                       G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
129{ 
130  //Build cross section table : Taken from Ref.7
131  G4double cs=0.0;
132  G4double kinEnergy = kineticEnergy;
133  if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
134  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
135
136  G4double value0,value1;
137  CalculateIntegrals(p,Z,kinEnergy,value0,value1);
138 
139  if(value1 > 0.0) cs = 1./value1;
140
141  return cs;
142} 
143//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144
145void 
146G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
147                                               G4double safety)
148{   
149  G4double kineticEnergy = dynParticle->GetKineticEnergy();
150  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)) return ;
151
152  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
153  G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
154  G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2;
155
156  ///////////////////////////////////////////
157  // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4
158  G4double  eloss = theManager->GetEnergy(particle,tPathLength,currentCouple);
159  if(eloss>0.5*kineticEnergy)eloss=kineticEnergy-eloss;//exchange possibility between target atomic e- and incident particle
160  G4double ee       = kineticEnergy - 0.5*eloss;
161  G4double ttau     = ee/electron_mass_c2;
162  G4double ttau2    = ttau*ttau;
163  G4double epsilonpp= eloss/ee;
164  G4double cst1=epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
165
166  kineticEnergy *= (1 - cst1);
167  ///////////////////////////////////////////
168  // additivity rule for mixture and compound xsection calculation
169  const G4Material* mat = currentCouple->GetMaterial();
170  G4int nelm = mat->GetNumberOfElements();
171  const G4ElementVector* theElementVector = mat->GetElementVector();
172  const G4double* theFraction = mat->GetFractionVector();
173  G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume();
174  G4double llambda0=0.,llambda1=0.;
175  for(G4int i=0;i<nelm;i++)
176    {
177      G4double l0,l1;
178      CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,l0,l1);
179      llambda0 += (theFraction[i]/l0);
180      llambda1 += (theFraction[i]/l1);
181    } 
182  if(llambda0>DBL_MIN) llambda0 =1./llambda0;
183  if(llambda1>DBL_MIN) llambda1 =1./llambda1;
184  G4double g1=0.0;
185  if(llambda1>DBL_MIN) g1 = llambda0/llambda1;
186
187  G4double x1,x0;
188  x0=g1/2.;
189  do
190    { 
191      x1 = x0-(x0*((1.+x0)*log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)
192      delta = std::abs( x1 - x0 );   
193      x0 = x1;  // new approximation becomes the old approximation for the next iteration
194    } while (delta > 1e-10);
195  G4double scrA = x1;
196
197  G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
198  G4double lambdan=0.;
199  G4bool mscatt=false,noscatt=false;
200
201  if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0;
202  if((lambdan<=1.0e-12))return;
203
204  G4double epsilon1=G4UniformRand();
205  G4double expn = exp(-lambdan);
206  if((epsilon1<expn)||insideskin)// no scattering
207    {noscatt=true;}
208  else if((epsilon1<((1.+lambdan)*expn)||(lambdan<1.)))
209    {
210      mscatt=false;
211      ws=G4UniformRand();
212      ws= 1.-2.*scrA*ws/(1.-ws + scrA);
213      if(acos(ws)<sqrt(scrA))//small angle approximation for theta less than screening angle
214      {G4int i=0;
215      do{i++;
216      ws=1.+0.5*atomPerVolume*tPathLength*log(G4UniformRand())/llambda1;
217      }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
218      if(i==19)ws=cos(sqrt(scrA));
219      }
220      G4double phi0=twopi*G4UniformRand(); 
221      us=sqrt(1.-ws*ws)*cos(phi0);
222      vs=sqrt(1.-ws*ws)*sin(phi0);
223      G4double rr=G4UniformRand();
224      x_coord=(rr*us);
225      y_coord=(rr*vs);
226      z_coord=((1.-rr)+rr*ws);
227    }
228  else
229    {
230      mscatt=true;
231      // Ref.2 subsection 4.4 "The best solution found"
232      // Sample first substep scattering angle
233      SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
234      phi1  = twopi*G4UniformRand();
235      cosPhi1 = cos(phi1);
236      sinPhi1 = sin(phi1);
237
238      // Sample second substep scattering angle
239      SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
240      phi2  = twopi*G4UniformRand();
241      cosPhi2 = cos(phi2);
242      sinPhi2 = sin(phi2);
243
244      // Scattering direction
245      us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
246      vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
247      ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
248    }
249   
250  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
251  G4ThreeVector newDirection(us,vs,ws);
252  newDirection.rotateUz(oldDirection);
253  fParticleChange->ProposeMomentumDirection(newDirection);
254 
255  if((safety > tlimitminfix)&&(latDisplasment))
256    { 
257
258      if(mscatt)
259        {
260          if(scrA<DBL_MIN)scrA=DBL_MIN;
261          if(llambda1<DBL_MIN)llambda1=DBL_MIN;
262          q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
263          if(q1<DBL_MIN)q1=DBL_MIN;
264          Gamma  = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
265          Eta    = atomPerVolume*tPathLength/llambda1;     
266          delta  = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
267
268          nu = G4UniformRand(); 
269          nu = std::sqrt(nu);
270          nu0 = (1.0 - nu)/2.;
271          nu1 = nu*delta;
272          nu2 = nu*(1.0-delta);
273          x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu0*us);
274          y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu0*vs);
275          z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu0*ws)  ;
276        }
277
278      // displacement is computed relatively to the end point
279      if(!noscatt)z_coord -= 1.0;//for noscatt zcoord z_coord !=0.
280      G4double rr = sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
281      G4double r  = rr*zPathLength;
282      /*
283      G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
284             << " sinTheta= " << sqrt(1.0 - ws*ws) << " r(mm)= " << r
285             << " trueStep(mm)= " << tPathLength
286             << " geomStep(mm)= " << zPathLength
287             << G4endl;
288      */
289      if(r > tlimitminfix) {
290
291        G4ThreeVector latDirection(x_coord/rr,y_coord/rr,z_coord/rr);
292        latDirection.rotateUz(oldDirection);
293
294        ComputeDisplacement(fParticleChange, latDirection, r, safety);
295      }     
296    }
297}
298
299//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
300
301void 
302G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
303                                                G4double &cost, G4double &sint)
304{
305  G4double u,Qn1,r1,tet;
306  G4double xi=0.;
307  Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
308
309if (Qn1<0.001)xi=-0.5*Qn1*log(G4UniformRand());
310else if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
311else
312{
313      // procedure described by Benedito in Ref.1
314      do{
315        r1=G4UniformRand();
316        u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
317        xi = 2.*u;
318        tet=acos(1.-xi);
319      }while(tet*r1*r1>sin(tet));
320}
321
322  if(xi<0.)xi=0.;
323  if(xi>2.)xi=2.;
324  cost=(1. - xi);
325  sint=sqrt(xi*(2.-xi));
326
327}
328
329//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
330// Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
331// linear log-log extrapolation between 1 GeV - 100 TeV
332
333void 
334G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 
335                                                 G4double kinEnergy,G4double &Lam0,
336                                                 G4double &Lam1)
337{ 
338  G4double summ00=0.0;
339  G4double summ10=0.0;
340  G4double x1,x2,y1,y2,acoeff,bcoeff;
341  G4double kineticE = kinEnergy;
342  if(kineticE<lowKEnergy)kineticE=lowKEnergy;
343  if(kineticE>highKEnergy)kineticE=highKEnergy;
344  kineticE /= eV;
345  G4double logE=log(kineticE);
346 
347  G4int  iZ = G4int(Z);
348  if(iZ > 103) iZ = 103;
349
350  G4int enerInd=0;
351  for(G4int i=1;i<106;i++)
352  {
353  if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
354  }
355
356  if(p==G4Electron::Electron())       
357    {
358    if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b
359      {
360        x1=ener[enerInd];
361        x2=ener[enerInd+1];       
362        y1=TCSE[iZ-1][enerInd];
363        y2=TCSE[iZ-1][enerInd+1];
364        acoeff=(y2-y1)/(x2*x2-x1*x1);
365        bcoeff=y2-acoeff*x2*x2;
366        summ00=acoeff*logE*logE+bcoeff;
367        summ00 =exp(summ00);
368        y1=FTCSE[iZ-1][enerInd];
369        y2=FTCSE[iZ-1][enerInd+1];
370        acoeff=(y2-y1)/(x2*x2-x1*x1);
371        bcoeff=y2-acoeff*x2*x2;
372        summ10=acoeff*logE*logE+bcoeff;
373        summ10 =exp(summ10);
374      }
375    else  //Interpolation of the form y=ax+b
376      { 
377        x1=ener[104];
378        x2=ener[105];       
379        y1=TCSE[iZ-1][104];
380        y2=TCSE[iZ-1][105];
381        summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
382        summ00 =exp(summ00);
383        y1=FTCSE[iZ-1][104];
384        y2=FTCSE[iZ-1][105];
385        summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
386        summ10 =exp(summ10);
387      }
388    }
389  if(p==G4Positron::Positron())       
390    {
391    if(kinEnergy<=1.0e+9)
392      {
393        x1=ener[enerInd];
394        x2=ener[enerInd+1];       
395        y1=TCSP[iZ-1][enerInd];
396        y2=TCSP[iZ-1][enerInd+1];
397        acoeff=(y2-y1)/(x2*x2-x1*x1);
398        bcoeff=y2-acoeff*x2*x2;
399        summ00=acoeff*logE*logE+bcoeff;
400        summ00 =exp(summ00);
401        y1=FTCSP[iZ-1][enerInd];
402        y2=FTCSP[iZ-1][enerInd+1];
403        acoeff=(y2-y1)/(x2*x2-x1*x1);
404        bcoeff=y2-acoeff*x2*x2;
405        summ10=acoeff*logE*logE+bcoeff;
406        summ10 =exp(summ10);
407      }
408    else
409      { 
410        x1=ener[104];
411        x2=ener[105];       
412        y1=TCSP[iZ-1][104];
413        y2=TCSP[iZ-1][105];
414        summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
415        summ00 =exp(summ00);
416        y1=FTCSP[iZ-1][104];
417        y2=FTCSP[iZ-1][105];
418        summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
419        summ10 =exp(summ10);
420      }
421    }
422   
423  summ00 *=barn;
424  summ10 *=barn;
425
426  Lam0=1./((1.+1./Z)*summ00);
427  Lam1=1./((1.+1./Z)*summ10);
428
429}
430
431//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
432//t->g->t step transformations taken from Ref.6
433
434G4double
435G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
436                                                         G4PhysicsTable* theTable,
437                                                         G4double currentMinimalStep)
438{
439  tPathLength = currentMinimalStep;
440  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
441  G4StepStatus stepStatus = sp->GetStepStatus();
442
443  const G4DynamicParticle* dp = track.GetDynamicParticle();
444
445  if(stepStatus == fUndefined) {
446    inside = false;
447    insideskin = false;
448    tlimit = geombig;
449    SetParticle( dp->GetDefinition() );
450  }
451
452  theLambdaTable = theTable;
453  currentCouple = track.GetMaterialCutsCouple();
454  currentMaterialIndex = currentCouple->GetIndex();
455  currentKinEnergy = dp->GetKineticEnergy();
456  currentRange = 
457    theManager->GetRangeFromRestricteDEDX(particle,currentKinEnergy,currentCouple);
458
459  lambda1 = GetLambda(currentKinEnergy);
460
461  // stop here if small range particle
462  if(inside) return tPathLength;           
463 
464  if(tPathLength > currentRange) tPathLength = currentRange;
465
466  G4double presafety = sp->GetSafety();
467
468  //G4cout << "G4GS::StepLimit tPathLength= "
469  //     <<tPathLength<<" safety= " << presafety
470  //       << " range= " <<currentRange<< " lambda= "<<lambda1
471  //     << " Alg: " << steppingAlgorithm <<G4endl;
472
473  // far from geometry boundary
474  if(currentRange < presafety)
475    {
476      inside = true;
477      return tPathLength; 
478    }
479
480  // standard  version
481  //
482  if (steppingAlgorithm == fUseDistanceToBoundary)
483    {
484      //compute geomlimit and presafety
485      G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
486   
487      // is far from boundary
488      if(currentRange <= presafety)
489        {
490          inside = true;
491          return tPathLength;   
492        }
493
494      smallstep += 1.;
495      insideskin = false;
496
497      if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
498        {
499          rangeinit = currentRange;
500          if(stepStatus == fUndefined) smallstep = 1.e10;
501          else  smallstep = 1.;
502
503          //define stepmin here (it depends on lambda!)
504          //rough estimation of lambda_elastic/lambda_transport
505          G4double rat = currentKinEnergy/MeV ;
506          rat = 1.e-3/(rat*(10.+rat)) ;
507          //stepmin ~ lambda_elastic
508          stepmin = rat*lambda1;
509          skindepth = skin*stepmin;
510          //define tlimitmin
511          tlimitmin = 10.*stepmin;
512          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
513
514          //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
515          //     << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
516          // constraint from the geometry
517          if((geomlimit < geombig) && (geomlimit > geommin))
518            {
519              if(stepStatus == fGeomBoundary) 
520                tgeom = geomlimit/facgeom;
521              else
522                tgeom = 2.*geomlimit/facgeom;
523            }
524            else
525              tgeom = geombig;
526
527        }
528
529      //step limit
530      tlimit = facrange*rangeinit;             
531      if(tlimit < facsafety*presafety)
532        tlimit = facsafety*presafety; 
533
534      //lower limit for tlimit
535      if(tlimit < tlimitmin) tlimit = tlimitmin;
536
537      if(tlimit > tgeom) tlimit = tgeom;
538
539      //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit 
540      //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
541
542      // shortcut
543      if((tPathLength < tlimit) && (tPathLength < presafety) &&
544         (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
545        return tPathLength;   
546
547      // step reduction near to boundary
548      if(smallstep < skin)
549        {
550          tlimit = stepmin;
551          insideskin = true;
552        }
553      else if(geomlimit < geombig)
554        {
555          if(geomlimit > skindepth)
556            {
557              if(tlimit > geomlimit-0.999*skindepth)
558                tlimit = geomlimit-0.999*skindepth;
559            }
560          else
561            {
562              insideskin = true;
563              if(tlimit > stepmin) tlimit = stepmin;
564            }
565        }
566
567      if(tlimit < stepmin) tlimit = stepmin;
568
569      if(tPathLength > tlimit) tPathLength = tlimit; 
570
571    }
572    // for 'normal' simulation with or without magnetic field
573    //  there no small step/single scattering at boundaries
574  else if(steppingAlgorithm == fUseSafety)
575    {
576      // compute presafety again if presafety <= 0 and no boundary
577      // i.e. when it is needed for optimization purposes
578      if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
579        presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
580
581      // is far from boundary
582      if(currentRange < presafety)
583        {
584          inside = true;
585          return tPathLength; 
586        }
587
588      if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
589        {
590          rangeinit = currentRange;
591          fr = facrange;
592          // 9.1 like stepping for e+/e- only (not for muons,hadrons)
593          if(mass < masslimite) 
594            {
595              if(lambda1 > currentRange)
596                rangeinit = lambda1;
597              if(lambda1 > lambdalimit)
598                fr *= 0.75+0.25*lambda1/lambdalimit;
599            }
600
601          //lower limit for tlimit
602          G4double rat = currentKinEnergy/MeV ;
603          rat = 1.e-3/(rat*(10.+rat)) ;
604          tlimitmin = 10.*lambda1*rat;
605          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
606        }
607      //step limit
608      tlimit = fr*rangeinit;               
609
610      if(tlimit < facsafety*presafety)
611        tlimit = facsafety*presafety;
612
613      //lower limit for tlimit
614      if(tlimit < tlimitmin) tlimit = tlimitmin;
615
616      if(tPathLength > tlimit) tPathLength = tlimit;
617    }
618 
619  // version similar to 7.1 (needed for some experiments)
620  else
621    {
622      if (stepStatus == fGeomBoundary)
623        {
624          if (currentRange > lambda1) tlimit = facrange*currentRange;
625          else                        tlimit = facrange*lambda1;
626
627          if(tlimit < tlimitmin) tlimit = tlimitmin;
628          if(tPathLength > tlimit) tPathLength = tlimit;
629        }
630    }
631  //G4cout << "tPathLength= " << tPathLength 
632  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
633  return tPathLength ;
634}
635
636//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
637
638G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
639{
640  par1 = -1. ; 
641  par2 = par3 = 0. ; 
642
643  //  do the true -> geom transformation
644  zPathLength = tPathLength;
645
646  // z = t for very small tPathLength
647  if(tPathLength < tlimitminfix) return zPathLength;
648
649  // this correction needed to run MSC with eIoni and eBrem inactivated
650  // and makes no harm for a normal run
651  if(tPathLength > currentRange)
652    tPathLength = currentRange ;
653
654  G4double tau   = tPathLength/lambda1 ;
655
656  if ((tau <= tausmall) || insideskin) {
657    zPathLength  = tPathLength;
658    if(zPathLength > lambda1) zPathLength = lambda1;
659    return zPathLength;
660  }
661
662  G4double zmean = tPathLength;
663  if (tPathLength < currentRange*dtrl) {
664    if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
665    else             zmean = lambda1*(1.-exp(-tau));
666  } else if(currentKinEnergy < mass) {
667    par1 = 1./currentRange ;
668    par2 = 1./(par1*lambda1) ;
669    par3 = 1.+par2 ;
670    if(tPathLength < currentRange)
671      zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
672    else
673      zmean = 1./(par1*par3) ;
674  } else {
675    G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,
676                                        currentCouple);
677
678    lambda11 = GetLambda(T1);
679
680    par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
681    par2 = 1./(par1*lambda1) ;
682    par3 = 1.+par2 ;
683    zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
684  }
685
686  zPathLength = zmean ;
687  //  sample z
688  if(samplez) {
689
690    const G4double  ztmax = 0.99, onethird = 1./3. ;
691    G4double zt = zmean/tPathLength ;
692
693    if (tPathLength > stepmin && zt < ztmax) {
694
695      G4double u,cz1;
696      if(zt >= onethird) {
697
698        G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
699        cz1 = 1.+cz ;
700        G4double u0 = cz/cz1 ;
701        G4double grej ;
702        do {
703          u = exp(log(G4UniformRand())/cz1) ;
704          grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
705        } while (grej < G4UniformRand()) ;
706
707      } else {
708        cz1 = 1./zt-1.;
709        u = 1.-exp(log(G4UniformRand())/cz1) ;
710      }
711      zPathLength = tPathLength*u ;
712    }
713  }
714  if(zPathLength > lambda1) zPathLength = lambda1;
715  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
716
717  return zPathLength;
718}
719
720//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
721
722G4double
723G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
724{
725  // step defined other than transportation
726  if(geomStepLength == zPathLength && tPathLength <= currentRange)
727    return tPathLength;
728
729  // t = z for very small step
730  zPathLength = geomStepLength;
731  tPathLength = geomStepLength;
732  if(geomStepLength < tlimitminfix) return tPathLength;
733 
734  // recalculation
735  if((geomStepLength > lambda1*tausmall) && !insideskin)
736    {
737      if(par1 <  0.)
738        tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
739      else 
740        {
741          if(par1*par3*geomStepLength < 1.)
742            tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
743          else 
744            tPathLength = currentRange;
745        } 
746    }
747  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
748  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
749
750  return tPathLength;
751}
752
753//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
754//Total & first transport x sections for e-/e+ generated from ELSEPA code
755
756void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
757{ 
758  G4String filename = "XSECTIONS.dat";
759
760  char* path = getenv("G4LEDATA");
761  if (!path)
762    {
763      G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set properly";
764      G4Exception(excep);
765    }
766
767  G4String pathString(path);
768  G4String dirFile = pathString + "/msc_GS/" + filename;
769  FILE *infile;
770  infile = fopen(dirFile,"r"); 
771  if (infile == 0)
772    {
773      G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
774      G4Exception(excep);
775    }
776
777  // Read parameters from tables and take logarithms
778  G4float aRead;
779  for(G4int i=0 ; i<106 ;i++){
780    fscanf(infile,"%f\t",&aRead);
781    if(aRead > 0.0) aRead = log(aRead);
782    else  aRead = 0.0;
783    ener[i]=aRead;
784  }       
785  for(G4int j=0;j<103;j++){
786    for(G4int i=0;i<106;i++){
787      fscanf(infile,"%f\t",&aRead);
788      if(aRead > 0.0) aRead = log(aRead);
789      else  aRead = 0.0;
790      TCSE[j][i]=aRead;
791    }       
792  }
793  for(G4int j=0;j<103;j++){
794    for(G4int i=0;i<106;i++){
795      fscanf(infile,"%f\t",&aRead);
796      if(aRead > 0.0) aRead = log(aRead);
797      else  aRead = 0.0;
798      FTCSE[j][i]=aRead;     
799    }       
800  }   
801  for(G4int j=0;j<103;j++){
802    for(G4int i=0;i<106;i++){
803      fscanf(infile,"%f\t",&aRead);
804      if(aRead > 0.0) aRead = log(aRead);
805      else  aRead = 0.0;
806      TCSP[j][i]=aRead;     
807    }       
808  }
809  for(G4int j=0;j<103;j++){
810    for(G4int i=0;i<106;i++){
811      fscanf(infile,"%f\t",&aRead);
812      if(aRead > 0.0) aRead = log(aRead);
813      else  aRead = 0.0;
814      FTCSP[j][i]=aRead;     
815    }       
816  }
817
818  fclose(infile);
819
820}
821
822//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.