G4GoudsmitSaundersonMscModel Class Reference

#include <G4GoudsmitSaundersonMscModel.hh>

Inheritance diagram for G4GoudsmitSaundersonMscModel:

G4VMscModel G4VEmModel

Public Member Functions

 G4GoudsmitSaundersonMscModel (const G4String &nam="GoudsmitSaunderson")
virtual ~G4GoudsmitSaundersonMscModel ()
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
void StartTracking (G4Track *)
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double, G4double, G4double)
virtual G4ThreeVector & SampleScattering (const G4ThreeVector &, G4double safety)
virtual G4double ComputeTruePathLengthLimit (const G4Track &track, G4double &currentMinimalStep)
virtual G4double ComputeGeomPathLength (G4double truePathLength)
virtual G4double ComputeTrueStepLength (G4double geomStepLength)

Detailed Description

Definition at line 73 of file G4GoudsmitSaundersonMscModel.hh.

Constructor & Destructor Documentation

G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel ( const G4String &  nam = "GoudsmitSaunderson"  ) 

Definition at line 97 of file G4GoudsmitSaundersonMscModel.cc.

References G4cout, G4endl, G4LossTableManager::Instance(), and G4VMscModel::samplez.

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 }

G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel (  )  [virtual]

Definition at line 123 of file G4GoudsmitSaundersonMscModel.cc.

00124 {
00125   delete GSTable;
00126 }

Member Function Documentation

G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition *  particle,
G4double  KineticEnergy,
G4double  AtomicNumber,
G4double  ,
G4double  ,
G4double   
) [virtual]

Reimplemented from G4VEmModel.

Definition at line 139 of file G4GoudsmitSaundersonMscModel.cc.

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 }

G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength ( G4double  truePathLength  )  [virtual]

Reimplemented from G4VMscModel.

Definition at line 647 of file G4GoudsmitSaundersonMscModel.cc.

References G4VMscModel::dtrl, G4UniformRand, G4VMscModel::GetEnergy(), G4VMscModel::GetTransportMeanFreePath(), and G4VMscModel::samplez.

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 }

G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit ( const G4Track &  track,
G4double &  currentMinimalStep 
) [virtual]

Reimplemented from G4VMscModel.

Definition at line 454 of file G4GoudsmitSaundersonMscModel.cc.

References G4VMscModel::ComputeGeomLimit(), G4VMscModel::ComputeSafety(), G4VMscModel::ConvertTrueToGeom(), G4VMscModel::facgeom, G4VMscModel::facrange, G4VMscModel::facsafety, fGeomBoundary, fUseDistanceToBoundary, fUseSafety, G4Track::GetDynamicParticle(), G4MaterialCutsCouple::GetIndex(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetMaterialCutsCouple(), G4Step::GetPreStepPoint(), G4VMscModel::GetRange(), G4Track::GetStep(), G4VMscModel::GetTransportMeanFreePath(), G4VEmModel::SetCurrentCouple(), G4VMscModel::skin, G4InuclParticleNames::sp, and G4VMscModel::steppingAlgorithm.

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 }

G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength ( G4double  geomStepLength  )  [virtual]

Reimplemented from G4VMscModel.

Definition at line 732 of file G4GoudsmitSaundersonMscModel.cc.

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 }

void G4GoudsmitSaundersonMscModel::Initialise ( const G4ParticleDefinition *  ,
const G4DataVector &   
) [virtual]

Implements G4VEmModel.

Definition at line 128 of file G4GoudsmitSaundersonMscModel.cc.

References G4VMscModel::GetParticleChangeForMSC(), and G4VMscModel::skin.

00130 { 
00131   skindepth=skin*stepmin;
00132   SetParticle(p);
00133   fParticleChange = GetParticleChangeForMSC(p);
00134 }

G4ThreeVector & G4GoudsmitSaundersonMscModel::SampleScattering ( const G4ThreeVector &  ,
G4double  safety 
) [virtual]

Reimplemented from G4VMscModel.

Definition at line 154 of file G4GoudsmitSaundersonMscModel.cc.

References G4VMscModel::dtrl, G4VMscModel::fDisplacement, G4UniformRand, G4VMscModel::GetDEDX(), G4Material::GetElementVector(), G4VMscModel::GetEnergy(), G4MaterialCutsCouple::GetMaterial(), G4Material::GetNumberOfElements(), G4Material::GetVecNbOfAtomsPerVolume(), G4ParticleChangeForMSC::ProposeMomentumDirection(), and G4InuclParticleNames::s0.

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 }

void G4GoudsmitSaundersonMscModel::StartTracking ( G4Track *   )  [virtual]

Reimplemented from G4VEmModel.

Definition at line 441 of file G4GoudsmitSaundersonMscModel.cc.

References G4DynamicParticle::GetDefinition(), and G4Track::GetDynamicParticle().

00442 {
00443   SetParticle(track->GetDynamicParticle()->GetDefinition());
00444   firstStep = true; 
00445   inside = false;
00446   insideskin = false;
00447   tlimit = geombig;
00448 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:06 2013 for Geant4 by  doxygen 1.4.7