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

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

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