G4GoudsmitSaundersonMscModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4GoudsmitSaundersonMscModel.cc 66592 2012-12-23 09:34:55Z vnivanch $
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 // File name:     G4GoudsmitSaundersonMscModel
00033 //
00034 // Author:        Omrane Kadri
00035 //
00036 // Creation date: 20.02.2009
00037 //
00038 // Modifications:
00039 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
00040 //
00041 // 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta 
00042 //                     sampling from SampleCosineTheta() which means the splitting 
00043 //                     step into two sub-steps occur only for msc regime
00044 //
00045 // 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
00046 //                     adding a theta min limit due to screening effect of the atomic nucleus
00047 // 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
00048 //                     within CalculateIntegrals method
00049 // 05.10.2009 O.Kadri: tuning small angle theta distributions
00050 //                     assuming the case of lambdan<1 as single scattering regime
00051 //                     tuning theta sampling for theta below the screening angle
00052 // 08.02.2010 O.Kadri: bugfix in compound xsection calculation and small angle computation
00053 //                     adding a rejection condition to hard collision angular sampling
00054 //                     ComputeTruePathLengthLimit was taken from G4WentzelVIModel
00055 // 26.03.2010 O.Kadri: direct xsection calculation not inverse of the inverse
00056 //                     angular sampling without large angle rejection method
00057 //                     longitudinal displacement is computed exactly from <z>
00058 // 12.05.2010 O.Kadri: exchange between target and projectile has as a condition the particle type (e-/e-)
00059 //                     some cleanup to minimize time consuming (adding lamdan12 & Qn12, changing the error to 1.0e-12 for scrA)
00060 //
00061 //
00062 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 
00063 //REFERENCES:
00064 //Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
00065 //Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
00066 //Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
00067 //Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
00068 //Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
00069 //Ref.6:G4UrbanMscModel G4 9.2; 
00070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00071 #include "G4GoudsmitSaundersonMscModel.hh"
00072 #include "G4GoudsmitSaundersonTable.hh"
00073 
00074 #include "G4PhysicalConstants.hh"
00075 #include "G4SystemOfUnits.hh"
00076 
00077 #include "G4ParticleChangeForMSC.hh"
00078 #include "G4MaterialCutsCouple.hh"
00079 #include "G4DynamicParticle.hh"
00080 #include "G4Electron.hh"
00081 #include "G4Positron.hh"
00082 
00083 #include "G4LossTableManager.hh"
00084 #include "G4Track.hh"
00085 #include "G4PhysicsTable.hh"
00086 #include "Randomize.hh"
00087 
00088 using namespace std;
00089 
00090 G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
00091 G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
00092 G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
00093 G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
00094 G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
00095 
00096 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00097 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
00098   : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
00099 { 
00100   currentKinEnergy=currentRange=skindepth=par1=par2=par3
00101     =zPathLength=truePathLength
00102     =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
00103     =lambda11=mass=0.0;
00104   currentMaterialIndex = -1;
00105 
00106   fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
00107   particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
00108   tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
00109   tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
00110   theManager=G4LossTableManager::Instance();
00111   inside=false;insideskin=false;
00112   samplez=false;
00113   firstStep = true; 
00114 
00115   GSTable = new G4GoudsmitSaundersonTable();
00116 
00117   if(ener[0] < 0.0){ 
00118     G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
00119     LoadELSEPAXSections();
00120   }
00121 }
00122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00123 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
00124 {
00125   delete GSTable;
00126 }
00127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00128 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
00129                                               const G4DataVector&)
00130 { 
00131   skindepth=skin*stepmin;
00132   SetParticle(p);
00133   fParticleChange = GetParticleChangeForMSC(p);
00134 }
00135 
00136 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00137 
00138 G4double 
00139 G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
00140                        G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
00141 {  
00142   G4double kinEnergy = kineticEnergy;
00143   if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
00144   if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
00145 
00146   G4double cs(0.0), cs0(0.0);
00147   CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
00148   
00149   return cs;
00150 }  
00151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00152 
00153 G4ThreeVector& 
00154 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection, G4double)
00155 {
00156   fDisplacement.set(0.0,0.0,0.0);
00157   G4double kineticEnergy = currentKinEnergy;
00158   //dynParticle->GetKineticEnergy();
00159   if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
00160      (tPathLength/tausmall < lambda1)) { return fDisplacement; }
00161 
00163   // Effective energy 
00164   G4double eloss  = 0.0;
00165   if (tPathLength > currentRange*dtrl) {
00166     eloss = kineticEnergy - 
00167       GetEnergy(particle,currentRange-tPathLength,currentCouple);
00168   } else {
00169     eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
00170   }
00171   /*
00172   G4double ttau      = kineticEnergy/electron_mass_c2;
00173   G4double ttau2     = ttau*ttau;
00174   G4double epsilonpp = eloss/kineticEnergy;
00175   G4double cst1  = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
00176   kineticEnergy *= (1 - cst1);
00177   */
00178   kineticEnergy -= 0.5*eloss;
00179  
00181   // additivity rule for mixture and compound xsection's
00182   const G4Material* mat = currentCouple->GetMaterial();
00183   const G4ElementVector* theElementVector = mat->GetElementVector();
00184   const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
00185   G4int nelm = mat->GetNumberOfElements();
00186   G4double s0(0.0), s1(0.0);
00187   lambda0 = 0.0;
00188   for(G4int i=0;i<nelm;i++)
00189     { 
00190       CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
00191       lambda0 += (theAtomNumDensityVector[i]*s0);
00192     } 
00193   if(lambda0>0.0) lambda0 =1./lambda0;
00194 
00195   // Newton-Raphson root's finding method of scrA from: 
00196   // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
00197   G4double g1=0.0;
00198   if(lambda1>0.0) { g1 = lambda0/lambda1; }
00199 
00200   G4double logx0,x1,delta;
00201   G4double x0=g1*0.5;
00202   // V.Ivanchenko added limit of the loop
00203   for(G4int i=0;i<1000;++i)
00204     {  
00205       logx0=std::log(1.+1./x0);
00206       x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
00207 
00208       // V.Ivanchenko cut step size of iterative procedure
00209       if(x1 < 0.0)         { x1 = 0.5*x0; }
00210       else if(x1 > 2*x0)   { x1 = 2*x0; }
00211       else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
00212       delta = std::fabs( x1 - x0 );    
00213       x0 = x1;
00214       if(delta < 1.0e-3*x1) { break;}
00215     }
00216   G4double scrA = x1;
00217 
00218   G4double lambdan=0.;
00219 
00220   if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
00221   if(lambdan<=1.0e-12) { return fDisplacement; }
00222  
00223   //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
00224   //    << " L1= " << lambda1 << G4endl;
00225 
00226   G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
00227   G4double Qn12 = 0.5*Qn1;
00228   
00229   G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
00230   G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
00231   G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
00232 
00233   G4double epsilon1=G4UniformRand();
00234   G4double expn = std::exp(-lambdan);
00235 
00236   if(epsilon1<expn)// no scattering 
00237     { return fDisplacement; }
00238   else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
00239     {
00240       G4double xi=G4UniformRand();
00241       xi= 2.*scrA*xi/(1.-xi + scrA);
00242       if(xi<0.)xi=0.;
00243       else if(xi>2.)xi=2.;      
00244       ws=(1. - xi);
00245       wss=std::sqrt(xi*(2.-xi));      
00246       G4double phi0=CLHEP::twopi*G4UniformRand(); 
00247       us=wss*cos(phi0);
00248       vs=wss*sin(phi0);
00249     }
00250   else // multiple scattering
00251     {
00252       // Ref.2 subsection 4.4 "The best solution found"
00253       // Sample first substep scattering angle
00254       SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
00255       G4double phi1  = CLHEP::twopi*G4UniformRand();
00256       cosPhi1 = cos(phi1);
00257       sinPhi1 = sin(phi1);
00258 
00259       // Sample second substep scattering angle
00260       SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
00261       G4double phi2  = CLHEP::twopi*G4UniformRand();
00262       cosPhi2 = cos(phi2);
00263       sinPhi2 = sin(phi2);
00264 
00265       // Overall scattering direction
00266       us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
00267       vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
00268       ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
00269       G4double sqrtA=sqrt(scrA);
00270       if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
00271       {
00272         G4int i=0;
00273         do{i++;
00274           ws=1.+Qn12*log(G4UniformRand());
00275         }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
00276         if(i>=19)ws=cos(sqrtA);
00277         wss=std::sqrt((1.-ws*ws));      
00278         us=wss*std::cos(phi1);
00279         vs=wss*std::sin(phi1);
00280       }
00281     }
00282     
00283   //G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
00284   G4ThreeVector newDirection(us,vs,ws);
00285   newDirection.rotateUz(oldDirection);
00286   fParticleChange->ProposeMomentumDirection(newDirection);
00287  
00288   // corresponding to error less than 1% in the exact formula of <z>
00289   if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
00290   else         { z_coord = (1.-std::exp(-Qn1))/Qn1; }
00291   G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
00292   x_coord  = rr*us;
00293   y_coord  = rr*vs;
00294 
00295   // displacement is computed relatively to the end point
00296   z_coord -= 1.0;
00297   z_coord *= zPathLength;
00298   /*
00299     G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
00300     << " sinTheta= " << sqrt(1.0 - ws*ws) 
00301     << " trueStep(mm)= " << tPathLength
00302     << " geomStep(mm)= " << zPathLength
00303     << G4endl;
00304   */
00305 
00306   fDisplacement.set(x_coord,y_coord,z_coord);
00307   fDisplacement.rotateUz(oldDirection);
00308 
00309   return fDisplacement;
00310 }     
00311 
00312 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00313 
00314 void 
00315 G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
00316                                                 G4double &cost, G4double &sint)
00317 {
00318   G4double r1,tet,xi=0.;
00319   G4double Qn1 = 2.* lambdan;
00320   if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*log(1.+1./scrA)-1.); }
00321   else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
00322   if (Qn1<0.001)
00323     {
00324       do{
00325         r1=G4UniformRand();
00326         xi=-0.5*Qn1*log(G4UniformRand());
00327         tet=acos(1.-xi);
00328       }while(tet*r1*r1>sin(tet));
00329     }
00330   else if(Qn1>0.5) { xi=2.*G4UniformRand(); }//isotropic distribution
00331   else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
00332 
00333 
00334   if(xi<0.)xi=0.;
00335   else if(xi>2.)xi=2.;
00336 
00337   cost=(1. - xi);
00338   sint=sqrt(xi*(2.-xi));
00339 
00340 }
00341 
00342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00343 // Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
00344 // linear log-log extrapolation between 1 GeV - 100 TeV
00345 
00346 void 
00347 G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 
00348                                                  G4double kinEnergy,G4double &Sig0,
00349                                                  G4double &Sig1)
00350 { 
00351   G4double x1,x2,y1,y2,acoeff,bcoeff;
00352   G4double kineticE = kinEnergy;
00353   if(kineticE<lowKEnergy)kineticE=lowKEnergy;
00354   if(kineticE>highKEnergy)kineticE=highKEnergy;
00355   kineticE /= eV;
00356   G4double logE=std::log(kineticE);
00357   
00358   G4int  iZ = G4int(Z);
00359   if(iZ > 103) iZ = 103;
00360 
00361   G4int enerInd=0;
00362   for(G4int i=0;i<105;i++)
00363   {
00364   if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
00365   }
00366 
00367   if(p==G4Electron::Electron())        
00368     {
00369     if(kineticE<=1.0e+9)//Interpolation of the form y=ax˛+b
00370       {
00371         x1=ener[enerInd];
00372         x2=ener[enerInd+1];       
00373         y1=TCSE[iZ-1][enerInd];
00374         y2=TCSE[iZ-1][enerInd+1];
00375         acoeff=(y2-y1)/(x2*x2-x1*x1);
00376         bcoeff=y2-acoeff*x2*x2;
00377         Sig0=acoeff*logE*logE+bcoeff;
00378         Sig0 =std::exp(Sig0);
00379         y1=FTCSE[iZ-1][enerInd];
00380         y2=FTCSE[iZ-1][enerInd+1];
00381         acoeff=(y2-y1)/(x2*x2-x1*x1);
00382         bcoeff=y2-acoeff*x2*x2;
00383         Sig1=acoeff*logE*logE+bcoeff;
00384         Sig1=std::exp(Sig1);
00385       }
00386     else  //Interpolation of the form y=ax+b
00387       {  
00388         x1=ener[104];
00389         x2=ener[105];       
00390         y1=TCSE[iZ-1][104];
00391         y2=TCSE[iZ-1][105];
00392         Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00393         Sig0=std::exp(Sig0);
00394         y1=FTCSE[iZ-1][104];
00395         y2=FTCSE[iZ-1][105];
00396         Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00397         Sig1=std::exp(Sig1);
00398       }
00399     }
00400   if(p==G4Positron::Positron())        
00401     {
00402     if(kinEnergy<=1.0e+9)
00403       {
00404         x1=ener[enerInd];
00405         x2=ener[enerInd+1];       
00406         y1=TCSP[iZ-1][enerInd];
00407         y2=TCSP[iZ-1][enerInd+1];
00408         acoeff=(y2-y1)/(x2*x2-x1*x1);
00409         bcoeff=y2-acoeff*x2*x2;
00410         Sig0=acoeff*logE*logE+bcoeff;
00411         Sig0 =std::exp(Sig0);
00412         y1=FTCSP[iZ-1][enerInd];
00413         y2=FTCSP[iZ-1][enerInd+1];
00414         acoeff=(y2-y1)/(x2*x2-x1*x1);
00415         bcoeff=y2-acoeff*x2*x2;
00416         Sig1=acoeff*logE*logE+bcoeff;
00417         Sig1=std::exp(Sig1);
00418       }
00419     else
00420       {  
00421         x1=ener[104];
00422         x2=ener[105];       
00423         y1=TCSP[iZ-1][104];
00424         y2=TCSP[iZ-1][105];
00425         Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00426         Sig0 =std::exp(Sig0);
00427         y1=FTCSP[iZ-1][104];
00428         y2=FTCSP[iZ-1][105];
00429         Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
00430         Sig1=std::exp(Sig1);
00431       }
00432     }
00433     
00434   Sig0 *= barn;
00435   Sig1 *= barn;
00436 
00437 }
00438 
00439 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00440 
00441 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track)
00442 {
00443   SetParticle(track->GetDynamicParticle()->GetDefinition());
00444   firstStep = true; 
00445   inside = false;
00446   insideskin = false;
00447   tlimit = geombig;
00448 }
00449 
00450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00451 //t->g->t step transformations taken from Ref.6 
00452 
00453 G4double 
00454 G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
00455                                                          G4double& currentMinimalStep)
00456 {
00457   tPathLength = currentMinimalStep;
00458   const G4DynamicParticle* dp = track.GetDynamicParticle();
00459   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
00460   G4StepStatus stepStatus = sp->GetStepStatus();
00461   currentCouple = track.GetMaterialCutsCouple();
00462   SetCurrentCouple(currentCouple); 
00463   currentMaterialIndex = currentCouple->GetIndex();
00464   currentKinEnergy = dp->GetKineticEnergy();
00465   currentRange = GetRange(particle,currentKinEnergy,currentCouple);
00466 
00467   lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
00468 
00469   // stop here if small range particle
00470   if(inside || tPathLength < tlimitminfix) {
00471     return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00472   }  
00473   if(tPathLength > currentRange) tPathLength = currentRange;
00474 
00475   G4double presafety = sp->GetSafety();
00476 
00477   //G4cout << "G4GS::StepLimit tPathLength= " 
00478   //     <<tPathLength<<" safety= " << presafety
00479   //       << " range= " <<currentRange<< " lambda= "<<lambda1
00480   //     << " Alg: " << steppingAlgorithm <<G4endl;
00481 
00482   // far from geometry boundary
00483   if(currentRange < presafety)
00484     {
00485       inside = true;
00486       return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
00487     }
00488 
00489   // standard  version
00490   //
00491   if (steppingAlgorithm == fUseDistanceToBoundary)
00492     {
00493       //compute geomlimit and presafety 
00494       G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
00495    
00496       // is far from boundary
00497       if(currentRange <= presafety)
00498         {
00499           inside = true;
00500           return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
00501         }
00502 
00503       smallstep += 1.;
00504       insideskin = false;
00505 
00506       if(firstStep || stepStatus == fGeomBoundary)
00507         {
00508           rangeinit = currentRange;
00509           if(firstStep) smallstep = 1.e10;
00510           else  smallstep = 1.;
00511 
00512           //define stepmin here (it depends on lambda!)
00513           //rough estimation of lambda_elastic/lambda_transport
00514           G4double rat = currentKinEnergy/MeV ;
00515           rat = 1.e-3/(rat*(10.+rat)) ;
00516           //stepmin ~ lambda_elastic
00517           stepmin = rat*lambda1;
00518           skindepth = skin*stepmin;
00519           //define tlimitmin
00520           tlimitmin = 10.*stepmin;
00521           if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00522 
00523           //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
00524           //     << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
00525           // constraint from the geometry 
00526           if((geomlimit < geombig) && (geomlimit > geommin))
00527             {
00528               if(stepStatus == fGeomBoundary)  
00529                 tgeom = geomlimit/facgeom;
00530               else
00531                 tgeom = 2.*geomlimit/facgeom;
00532             }
00533             else
00534               tgeom = geombig;
00535 
00536         }
00537 
00538       //step limit 
00539       tlimit = facrange*rangeinit;              
00540       if(tlimit < facsafety*presafety)
00541         tlimit = facsafety*presafety; 
00542 
00543       //lower limit for tlimit
00544       if(tlimit < tlimitmin) tlimit = tlimitmin;
00545 
00546       if(tlimit > tgeom) tlimit = tgeom;
00547 
00548       //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
00549       //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
00550 
00551       // shortcut
00552       if((tPathLength < tlimit) && (tPathLength < presafety) &&
00553          (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
00554         return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
00555 
00556       // step reduction near to boundary
00557       if(smallstep < skin)
00558         {
00559           tlimit = stepmin;
00560           insideskin = true;
00561         }
00562       else if(geomlimit < geombig)
00563         {
00564           if(geomlimit > skindepth)
00565             {
00566               if(tlimit > geomlimit-0.999*skindepth)
00567                 tlimit = geomlimit-0.999*skindepth;
00568             }
00569           else
00570             {
00571               insideskin = true;
00572               if(tlimit > stepmin) tlimit = stepmin;
00573             }
00574         }
00575 
00576       if(tlimit < stepmin) tlimit = stepmin;
00577 
00578       if(tPathLength > tlimit) tPathLength = tlimit; 
00579 
00580     }
00581     // for 'normal' simulation with or without magnetic field 
00582     //  there no small step/single scattering at boundaries
00583   else if(steppingAlgorithm == fUseSafety)
00584     {
00585       // compute presafety again if presafety <= 0 and no boundary
00586       // i.e. when it is needed for optimization purposes
00587       if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
00588         presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
00589 
00590       // is far from boundary
00591       if(currentRange < presafety)
00592         {
00593           inside = true;
00594           return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
00595         }
00596 
00597       if(firstStep || stepStatus == fGeomBoundary)
00598         {
00599           rangeinit = currentRange;
00600           fr = facrange;
00601           // 9.1 like stepping for e+/e- only (not for muons,hadrons)
00602           if(mass < masslimite) 
00603             {
00604               if(lambda1 > currentRange)
00605                 rangeinit = lambda1;
00606               if(lambda1 > lambdalimit)
00607                 fr *= 0.75+0.25*lambda1/lambdalimit;
00608             }
00609 
00610           //lower limit for tlimit
00611           G4double rat = currentKinEnergy/MeV ;
00612           rat = 1.e-3/(rat*(10.+rat)) ;
00613           tlimitmin = 10.*lambda1*rat;
00614           if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00615         }
00616       //step limit
00617       tlimit = fr*rangeinit;               
00618 
00619       if(tlimit < facsafety*presafety)
00620         tlimit = facsafety*presafety;
00621 
00622       //lower limit for tlimit
00623       if(tlimit < tlimitmin) tlimit = tlimitmin;
00624 
00625       if(tPathLength > tlimit) tPathLength = tlimit;
00626     }
00627   
00628   // version similar to 7.1 (needed for some experiments)
00629   else
00630     {
00631       if (stepStatus == fGeomBoundary)
00632         {
00633           if (currentRange > lambda1) tlimit = facrange*currentRange;
00634           else                        tlimit = facrange*lambda1;
00635 
00636           if(tlimit < tlimitmin) tlimit = tlimitmin;
00637           if(tPathLength > tlimit) tPathLength = tlimit;
00638         }
00639     }
00640   //G4cout << "tPathLength= " << tPathLength  
00641   // << " currentMinimalStep= " << currentMinimalStep << G4endl;
00642   return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00643 }
00644 
00645 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00646 // taken from Ref.6
00647 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
00648 {
00649   firstStep = false; 
00650   par1 = -1. ;  
00651   par2 = par3 = 0. ;  
00652 
00653   //  do the true -> geom transformation
00654   zPathLength = tPathLength;
00655 
00656   // z = t for very small tPathLength
00657   if(tPathLength < tlimitminfix) { return zPathLength; }
00658 
00659   // this correction needed to run MSC with eIoni and eBrem inactivated
00660   // and makes no harm for a normal run
00661   if(tPathLength > currentRange)
00662     { tPathLength = currentRange; }
00663 
00664   G4double tau   = tPathLength/lambda1 ;
00665 
00666   if ((tau <= tausmall) || insideskin) {
00667     zPathLength  = tPathLength;
00668     if(zPathLength > lambda1) { zPathLength = lambda1; }
00669     return zPathLength;
00670   }
00671 
00672   G4double zmean = tPathLength;
00673   if (tPathLength < currentRange*dtrl) {
00674     if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
00675     else             zmean = lambda1*(1.-exp(-tau));
00676   } else if(currentKinEnergy < mass || tPathLength == currentRange) {
00677     par1 = 1./currentRange ;
00678     par2 = 1./(par1*lambda1) ;
00679     par3 = 1.+par2 ;
00680     if(tPathLength < currentRange)
00681       zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
00682     else
00683       zmean = 1./(par1*par3) ;
00684   } else {
00685     G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
00686 
00687     lambda11 = GetTransportMeanFreePath(particle,T1);
00688 
00689     par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
00690     par2 = 1./(par1*lambda1) ;
00691     par3 = 1.+par2 ;
00692     zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
00693   }
00694 
00695   zPathLength = zmean ;
00696   //  sample z
00697   if(samplez) {
00698 
00699     const G4double  ztmax = 0.99;
00700     G4double zt = zmean/tPathLength ;
00701 
00702     if (tPathLength > stepmin && zt < ztmax) {
00703 
00704       G4double u,cz1;
00705       if(zt >= 0.333333333) {
00706 
00707         G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
00708         cz1 = 1.+cz ;
00709         G4double u0 = cz/cz1 ;
00710         G4double grej ;
00711         do {
00712           u = exp(log(G4UniformRand())/cz1) ;
00713           grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
00714         } while (grej < G4UniformRand()) ;
00715 
00716       } else {
00717         cz1 = 1./zt-1.;
00718         u = 1.-exp(log(G4UniformRand())/cz1) ;
00719       }
00720       zPathLength = tPathLength*u ;
00721     }
00722   }
00723   if(zPathLength > lambda1) zPathLength = lambda1;
00724   //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
00725 
00726   return zPathLength;
00727 }
00728 
00729 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00730 // taken from Ref.6
00731 G4double 
00732 G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
00733 {
00734   // step defined other than transportation 
00735   if(geomStepLength == zPathLength && tPathLength <= currentRange)
00736     return tPathLength;
00737 
00738   // t = z for very small step
00739   zPathLength = geomStepLength;
00740   tPathLength = geomStepLength;
00741   if(geomStepLength < tlimitminfix) return tPathLength;
00742   
00743   // recalculation
00744   if((geomStepLength > lambda1*tausmall) && !insideskin)
00745     {
00746       if(par1 <  0.)
00747         tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
00748       else 
00749         {
00750           if(par1*par3*geomStepLength < 1.)
00751             tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
00752           else 
00753             tPathLength = currentRange;
00754         }  
00755     }
00756   if(tPathLength < geomStepLength) tPathLength = geomStepLength;
00757   //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
00758 
00759   return tPathLength;
00760 }
00761 
00762 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00763 //Total & first transport x sections for e-/e+ generated from ELSEPA code
00764 
00765 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
00766 { 
00767   G4String filename = "XSECTIONS.dat";
00768 
00769   char* path = getenv("G4LEDATA");
00770   if (!path)
00771     {
00772       G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
00773                   FatalException,
00774                   "Environment variable G4LEDATA not defined");
00775       return;
00776     }
00777 
00778   G4String pathString(path);
00779   G4String dirFile = pathString + "/msc_GS/" + filename;
00780   FILE *infile;
00781   infile = fopen(dirFile,"r"); 
00782   if (infile == 0)
00783     {
00784       G4ExceptionDescription ed;
00785       ed << "Data file <" + dirFile + "> is not opened!" << G4endl;
00786       G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00787                   "em0003",FatalException,ed);
00788       return;
00789     }
00790 
00791   // Read parameters from tables and take logarithms
00792   G4float aRead;
00793   for(G4int i=0 ; i<106 ;i++){
00794     if(1 == fscanf(infile,"%f\t",&aRead)) {
00795       if(aRead > 0.0) { aRead = log(aRead); }
00796       else            { aRead = 0.0; }
00797     } else {
00798       G4ExceptionDescription ed;
00799       ed << "Error reading <" + dirFile + "> loop #1 i= " << i << G4endl;
00800       G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00801                   "em0003",FatalException,ed);
00802       return;
00803     }
00804     ener[i]=aRead;
00805   }        
00806   for(G4int j=0;j<103;j++){
00807     for(G4int i=0;i<106;i++){
00808       if(1 == fscanf(infile,"%f\t",&aRead)) {
00809         if(aRead > 0.0) { aRead = log(aRead); }
00810         else            { aRead = 0.0; }
00811       } else {
00812         G4ExceptionDescription ed;
00813         ed << "Error reading <" + dirFile + "> loop #2 j= " << j 
00814            << "; i= " << i << G4endl;
00815         G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00816                     "em0003",FatalException,ed);
00817         return;
00818       }
00819       TCSE[j][i]=aRead;
00820     }        
00821   }
00822   for(G4int j=0;j<103;j++){
00823     for(G4int i=0;i<106;i++){
00824       if(1 == fscanf(infile,"%f\t",&aRead)) {
00825         if(aRead > 0.0) { aRead = log(aRead); }
00826         else            { aRead = 0.0; }
00827       } else {
00828         G4ExceptionDescription ed;
00829         ed << "Error reading <" + dirFile + "> loop #3 j= " << j 
00830            << "; i= " << i << G4endl;
00831         G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00832                     "em0003",FatalException,ed);
00833         return;
00834       }
00835       FTCSE[j][i]=aRead;      
00836     }        
00837   }    
00838   for(G4int j=0;j<103;j++){
00839     for(G4int i=0;i<106;i++){
00840       if(1 == fscanf(infile,"%f\t",&aRead)) {
00841         if(aRead > 0.0) { aRead = log(aRead); }
00842         else            { aRead = 0.0; }
00843       } else {
00844         G4ExceptionDescription ed;
00845         ed << "Error reading <" + dirFile + "> loop #4 j= " << j 
00846            << "; i= " << i << G4endl;
00847         G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00848                     "em0003",FatalException,ed);
00849         return;
00850       }
00851       TCSP[j][i]=aRead;      
00852     }        
00853   }
00854   for(G4int j=0;j<103;j++){
00855     for(G4int i=0;i<106;i++){
00856       if(1 == fscanf(infile,"%f\t",&aRead)) {
00857         if(aRead > 0.0) { aRead = log(aRead); }
00858         else            { aRead = 0.0; }
00859       } else {
00860         G4ExceptionDescription ed;
00861         ed << "Error reading <" + dirFile + "> loop #5 j= " << j 
00862            << "; i= " << i << G4endl;
00863         G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
00864                     "em0003",FatalException,ed);
00865         return;
00866       }
00867       FTCSP[j][i]=aRead;      
00868     }        
00869   }
00870 
00871   fclose(infile);
00872 
00873 }
00874 
00875 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:48:24 2013 for Geant4 by  doxygen 1.4.7