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

Last change on this file since 1315 was 1315, checked in by garnier, 14 years ago

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