Index: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WentzelVIModel.cc,v 1.16 2008/11/19 11:47:50 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4WentzelVIModel.cc,v 1.32 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -52,5 +52,4 @@
 //
 
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -60,6 +59,4 @@
 #include "G4LossTableManager.hh"
 #include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
 #include "G4PhysicsTableHelper.hh"
 #include "G4ElementVector.hh"
@@ -71,4 +68,7 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
+G4double G4WentzelVIModel::FormFactor[]    = {0.0};
 
 using namespace std;
@@ -96,8 +96,7 @@
   thePositron = G4Positron::Positron();
   theProton   = G4Proton::Proton();
-  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
+  lowEnergyLimit = 0.1*keV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff  = twopi*p0*p0;
-  constn = 6.937e-6/(MeV*MeV);
   tkin = targetZ = mom2 = DBL_MIN;
   ecut = etag = DBL_MAX;
@@ -106,7 +105,20 @@
   xsecn.resize(nelments);
   prob.resize(nelments);
-  for(size_t j=0; j<100; j++) {
-    FF[j]    = 0.0;
-  } 
+
+  // Thomas-Fermi screening radii
+  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
+
+  if(0.0 == ScreenRSquare[0]) {
+    G4double a0 = electron_mass_c2/0.88534; 
+    G4double constn = 6.937e-6/(MeV*MeV);
+
+    ScreenRSquare[0] = alpha2*a0*a0;
+    for(G4int j=1; j<100; j++) {
+      G4double x = a0*fNistManager->GetZ13(j);
+      ScreenRSquare[j] = alpha2*x*x;
+      x = fNistManager->GetA27(j); 
+      FormFactor[j] = constn*x*x;
+    } 
+  }
 }
 
@@ -132,13 +144,6 @@
   if(!isInitialized) {
     isInitialized = true;
-
-    if (pParticleChange)
-      fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-    else
-      fParticleChange = new G4ParticleChangeForMSC();
-
-    safetyHelper = G4TransportationManager::GetTransportationManager()
-      ->GetSafetyHelper();
-    safetyHelper->InitialiseHelper();
+    fParticleChange = GetParticleChangeForMSC();
+    InitialiseSafetyHelper();
   }
 }
@@ -156,5 +161,5 @@
   SetupKinematic(ekin, cutEnergy);
   SetupTarget(Z, ekin);
-  G4double xsec = ComputeTransportXSectionPerVolume();
+  G4double xsec = ComputeTransportXSectionPerAtom();
   /*  
   G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2 
@@ -167,5 +172,5 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume()
+G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
 {
   G4double xSection = 0.0;
@@ -189,26 +194,27 @@
       y = 0.0;
     }
-    xSection += y/targetZ;
-  }
-  /*
-  G4cout << "G4WentzelVIModel:XS per A " << " Z= " << Z << " e(MeV)= " << kinEnergy/MeV 
+    xSection = y;
+  }
+  /*  
+  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
+	 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
 	 << " cut(MeV)= " << ecut/MeV  
   	 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
-	 << " zmaxN= " << (1.0 - cosTetMsxNuc)/screenZ << G4endl;
+	 << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ 
+         << " costm= " << cosTetMaxNuc2 << G4endl;
   */
-
   // scattering off nucleus
   if(cosTetMaxNuc2 < 1.0) {
     x  = 1.0 - cosTetMaxNuc2;
     x1 = screenZ*formfactA;
-    x2 = 1.0/(1.0 - x1); 
+    x2 = 1.0 - x1; 
     x3 = x/screenZ;
     x4 = formfactA*x;
     // low-energy limit
     if(x3 < numlimit && x1 < numlimit) {
-      y = 0.5*x3*x3*x2*x2*x2*(1.0 - 1.333333*x3 + 1.5*x3*x3 - 1.5*x1
-      			      + 3.0*x1*x1 + 2.666666*x3*x1);
+      y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1
+		     + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2);
       // high energy limit
-    } else if(1.0 < x1) {
+    } else if(x2 <= 0.0) {
       x4 = x1*(1.0 + x3);
       y  = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
@@ -216,6 +222,8 @@
     } else {
       y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4)) 
-	   - x3/(1. + x3) - x4/(1. + x4))*x2*x2; 
-    }
+	   - x3/(1. + x3) - x4/(1. + x4))/(x2*x2); 
+    }
+    //G4cout << "y= " << y << " x1= " <<x1<<"  x2= " <<x2
+    //	   <<"  x3= "<<x3<<"  x4= " << x4<<G4endl;
     if(y < 0.0) {
       nwarnings++;
@@ -230,9 +238,15 @@
       y = 0.0;
     }
-    xSection += y; 
-  }
-  xSection *= (coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2); 
-  //  G4cout << "   XStotal= " << xSection/barn << " screenZ= " << screenZ 
-  //	 << " formF= " << formfactA << " for " << p->GetParticleName() << G4endl;
+    xSection += y*targetZ; 
+  }
+  xSection *= kinFactor;
+  /*
+  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
+	 << " screenZ= " << screenZ << " formF= " << formfactA 
+	 << " for " << particle->GetParticleName() 
+  	 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
+	 << " x= " << x 
+	 << G4endl;
+  */
   return xSection; 
 }
@@ -277,5 +291,5 @@
   // i.e. when it is needed for optimization purposes
   if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 
-    presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
+    presafety = ComputeSafety(sp->GetPosition(), tlimit); 
   /*
   G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= " 
@@ -291,5 +305,6 @@
     G4double rlimit = facrange*lambda0;
     G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
-    if(rcut > rlimit) rlimit = std::pow(2.0*rcut*rcut*lambda0,0.33333333);
+    if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
+    rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
     if(rlimit < tlimit) tlimit = rlimit;
   }
@@ -404,5 +419,5 @@
   } else {
 
-    // define threshold angle as 2 sigma of central value
+    // define threshold angle between single and multiple scattering 
     cosThetaMin = 1.0 - 3.0*x1;
 
@@ -570,35 +585,6 @@
 
     if(r > tlimitminfix) {
-      G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-      G4double fac= 1.;
-      if(r >= safety) {
-	//  ******* so safety is computed at boundary too ************
-	G4double newsafety = 
-	  safetyHelper->ComputeSafety(Position) - tlimitminfix;
-        if(newsafety <= 0.0) fac = 0.0;
-	else if(r > newsafety) fac = newsafety/r ;
-        //G4cout << "NewSafety= " << newsafety << " fac= " << fac 
-	// << " r= " << r << " sint= " << sint << " pos " << Position << G4endl;
-      }  
-
-      if(fac > 0.) {
-	// compute new endpoint of the Step
-	G4ThreeVector newPosition = Position + fac*pos;
-
-	// check safety after displacement
-	G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	// displacement to boundary
-	if(postsafety <= 0.0) {
-	  safetyHelper->Locate(newPosition, newDirection);
-
-	  // not on the boundary
-	} else { 
-	  safetyHelper->ReLocateWithinVolume(newPosition);
-	  // if(fac < 1.0) G4cout << "NewPosition " << newPosition << G4endl;
-	}
-     
-	fParticleChange->ProposePosition(newPosition);
-      } 
+      pos /= r;
+      ComputeDisplacement(fParticleChange, pos, r, safety);
     }
   }
@@ -624,6 +610,4 @@
   G4double xs = 0.0;
 
-  G4double fac = coeff*chargeSquare*invbeta2/mom2;
-
   for (G4int i=0; i<nelm; i++) {
     SetupTarget((*theElementVector)[i]->GetZ(), tkin);
@@ -635,5 +619,5 @@
     cosTetMaxNuc2  = std::max(cosnm,cosThetaMin); 
     cosTetMaxElec2 = std::max(cosem,cosThetaMin); 
-    xtsec += ComputeTransportXSectionPerVolume()*density;
+    xtsec += ComputeTransportXSectionPerAtom()*density;
     // return limit back
     cosTetMaxElec2 = cosem;
@@ -643,5 +627,5 @@
     G4double nsec = 0.0;
     G4double x1 = 1.0 - cosThetaMin + screenZ;
-    G4double f  = fac*targetZ*density; 
+    G4double f  = kinFactor*density; 
 
     // scattering off electrons
@@ -656,5 +640,6 @@
       G4double s  = screenZ*formfactA;
       G4double z1 = 1.0 - cosnm + screenZ;
-      G4double d  = (1.0 - s)/formfactA;
+      G4double s1 = 1.0 - s;
+      G4double d  = s1/formfactA;
 
       // check numerical limit
@@ -667,6 +652,5 @@
 	G4double x2 = x1 + d;
 	G4double z2 = z1 + d;
-	nsec = (1.0 + 2.0*s)*((cosThetaMin - cosnm)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
-			   2.0*log(z1*x2/(z2*x1))/d);
+	nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
       }
       nsec *= f*targetZ;
@@ -801,5 +785,6 @@
       G4double mom22 = t1*(t1 + 2.0*mass);
       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
-      if(ctm < 1.0) cosTetMaxElec = ctm;
+      if(ctm <  1.0) cosTetMaxElec = ctm;
+      if(ctm < -1.0) cosTetMaxElec = -1.0;
     }
   }
@@ -808,11 +793,2 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4WentzelVIModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
-                                     const G4MaterialCutsCouple*,
-				     const G4DynamicParticle*,
-				     G4double,
-				     G4double)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
