Index: trunk/source/processes/electromagnetic/standard/src/G4ASTARStopping.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ASTARStopping.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4ASTARStopping.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4ASTARStopping.cc,v 1.8 2008/11/24 18:28:09 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 
 //---------------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BetheBlochModel.cc,v 1.31 2009/04/24 14:24:00 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4BetheBlochModel.cc,v 1.34 2009/11/11 23:22:27 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -234,9 +234,9 @@
   G4double eexc  = material->GetIonisation()->GetMeanExcitationEnergy();
   G4double eexc2 = eexc*eexc;
-  G4double cden  = material->GetIonisation()->GetCdensity();
-  G4double mden  = material->GetIonisation()->GetMdensity();
-  G4double aden  = material->GetIonisation()->GetAdensity();
-  G4double x0den = material->GetIonisation()->GetX0density();
-  G4double x1den = material->GetIonisation()->GetX1density();
+  //G4double cden  = material->GetIonisation()->GetCdensity();
+  //G4double mden  = material->GetIonisation()->GetMdensity();
+  //G4double aden  = material->GetIonisation()->GetAdensity();
+  //G4double x0den = material->GetIonisation()->GetX0density();
+  //G4double x1den = material->GetIonisation()->GetX1density();
 
   G4double eDensity = material->GetElectronDensity();
@@ -252,8 +252,9 @@
   // density correction
   G4double x = log(bg2)/twoln10;
-  if ( x >= x0den ) {
-    dedx -= twoln10*x - cden ;
-    if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
-  }
+  //if ( x >= x0den ) {
+  //  dedx -= twoln10*x - cden ;
+  //  if ( x < x1den ) dedx -= aden*pow((x1den-x),mden) ;
+  //}
+  dedx -= material->GetIonisation()->DensityCorrection(x);
 
   // shell correction
@@ -283,11 +284,11 @@
 					     G4double length)
 {
-  const G4ParticleDefinition* p = dp->GetDefinition();
-  const G4Material* mat = couple->GetMaterial();
-  G4double preKinEnergy = dp->GetKineticEnergy();
-  G4double e = preKinEnergy - eloss*0.5;
-  if(e < 0.0) e = preKinEnergy*0.5;
-
   if(isIon) {
+    const G4ParticleDefinition* p = dp->GetDefinition();
+    const G4Material* mat = couple->GetMaterial();
+    G4double preKinEnergy = dp->GetKineticEnergy();
+    G4double e = preKinEnergy - eloss*0.5;
+    if(e < 0.0) e = preKinEnergy*0.5;
+
     G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
     GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
@@ -295,5 +296,5 @@
     eloss += length*corr->IonHighOrderCorrections(p,couple,e);
   }
-
+  /*
   if(nuclearStopping && preKinEnergy*proton_mass_c2/mass < chargeSquare*100.*MeV) {
 
@@ -307,12 +308,12 @@
       eloss += nloss;
     }
-    /*
+    
     G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
     	   << " de= " << eloss << " NIEL= " << nloss 
 	   << " dynQ= " << dp->GetCharge()/eplus << G4endl;
-    */
+    
     fParticleChange->ProposeNonIonizingEnergyDeposit(nloss);
   }
-
+  */
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4BetheHeitlerModel.cc,v 1.13 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BohrFluctuations.cc,v 1.7 2009/02/19 19:17:50 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4BohrFluctuations.cc,v 1.8 2009/09/29 11:33:22 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -63,5 +63,5 @@
  :G4VEmFluctuationModel(nam),
   particle(0),
-  minNumberInteractionsBohr(10.0),
+  minNumberInteractionsBohr(2.0),
   minFraction(0.2),
   xmin(0.2),
@@ -110,6 +110,6 @@
     }
     siga = sqrt(siga);
-    //G4cout << "siga= " << siga << G4endl;
     G4double twomeanLoss = meanLoss + meanLoss;
+    //G4cout << "siga= " << siga << " 2edp= " << twomeanLoss <<G4endl;
 
     if(twomeanLoss < siga) {
@@ -130,5 +130,5 @@
     loss = meanLoss*n/navr;
   }
-  //  G4cout << "loss= " << loss << G4endl;
+  //G4cout << "loss= " << loss << G4endl;
 
   return loss;
Index: trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggIonModel.cc,v 1.24 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4BraggIonModel.cc,v 1.26 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -136,4 +136,5 @@
 					       G4double kineticEnergy)
 {
+  //G4cout << "G4BraggIonModel::GetChargeSquareRatio e= " <<  kineticEnergy << G4endl;
   // this method is called only for ions
   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,kineticEnergy);
@@ -148,4 +149,5 @@
 					    G4double kineticEnergy)
 {
+  //G4cout << "G4BraggIonModel::GetParticleCharge e= " <<  kineticEnergy << G4endl;
   // this method is called only for ions
   return corr->GetParticleCharge(p,mat,kineticEnergy);
@@ -252,5 +254,5 @@
 					   G4double& eloss,
 					   G4double&,
-					   G4double length)
+					   G4double /*length*/)
 {
   // this method is called only for ions
@@ -261,8 +263,10 @@
   if(e < 0.0) e = preKinEnergy*0.5;
 
+  //G4cout << "G4BraggIonModel::CorrectionsAlongStep e= " <<  e << G4endl;
+
   G4double q2 = corr->EffectiveChargeSquareRatio(p,mat,e);
   GetModelOfFluctuations()->SetParticleAndCharge(p, q2);
   eloss *= q2*corr->EffectiveChargeCorrection(p,mat,e)/corrFactor; 
-
+  /*
   if(nuclearStopping) {
 
@@ -276,11 +280,13 @@
       eloss += nloss;
     }
-    /*
+    
     G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
     	   << " de= " << eloss << " NIEL= " << nloss 
 	   << " dynQ= " << dp->GetCharge()/eplus << G4endl;
-    */
+    
     fParticleChange->ProposeNonIonizingEnergyDeposit(nloss);
   }
+  */
+  //G4cout << "G4BraggIonModel::CorrectionsAlongStep end" << G4endl;
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4BraggModel.cc,v 1.23 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -240,5 +240,5 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
+/*
 void G4BraggModel::CorrectionsAlongStep(const G4MaterialCutsCouple* couple,
 					const G4DynamicParticle* dp,
@@ -263,12 +263,13 @@
       eloss += nloss;
     }
-    /*
+    
     G4cout << "G4ionIonisation::CorrectionsAlongStep: e= " << preKinEnergy
     	   << " de= " << eloss << " NIEL= " << nloss 
 	   << " dynQ= " << dp->GetCharge()/eplus << G4endl;
-    */
+    
     fParticleChange->ProposeNonIonizingEnergyDeposit(nloss);
   }
 }
+*/
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4ComptonScattering.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
Index: trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering52.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering52.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering52.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4ComptonScattering52.cc,v 1.7 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
Index: trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScattering.cc,v 1.20 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4CoulombScattering.cc,v 1.25 2009/10/28 10:14:13 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -52,5 +52,7 @@
 #include "G4CoulombScatteringModel.hh"
 #include "G4eCoulombScatteringModel.hh"
+//#include "G4hCoulombScatteringModel.hh"
 #include "G4Electron.hh"
+#include "G4Proton.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -103,9 +105,9 @@
     aParticle = p;
     G4double mass = p->GetPDGMass();
+    G4String name = p->GetParticleName();
     if (mass > GeV || p->GetParticleType() == "nucleus") {
       SetBuildTableFlag(false);
       verboseLevel = 0;
     } else {
-      G4String name = p->GetParticleName();
       if(name != "e-" && name != "e+" &&
          name != "mu+" && name != "mu-" && name != "pi+" && 
@@ -115,4 +117,12 @@
     G4double emin = MinKinEnergy();
     G4double emax = MaxKinEnergy();
+    //    G4gCoulombScatteringModel* model = new G4hCoulombScatteringModel();
+    G4eCoulombScatteringModel* model = new G4eCoulombScatteringModel();
+    model->SetPolarAngleLimit(PolarAngleLimit());
+    model->SetLowEnergyLimit(emin);
+    model->SetHighEnergyLimit(emax);
+    AddEmModel(1, model);
+    /*
+
     G4double eth  = thEnergy;
     if(mass < MeV) eth  = thEnergyElec; 
@@ -131,4 +141,5 @@
       AddEmModel(2, model);
     }
+    */
   }
 }
@@ -139,6 +150,7 @@
 {
   G4cout << "      " << PolarAngleLimit()/degree
-	 << " < Theta(degree) < 180" 
+	 << " < Theta(degree) < 180"
 	 << ", Eth(MeV)= ";
+
   if(aParticle->GetPDGMass() < MeV) G4cout << thEnergyElec;
   else                              G4cout << thEnergy;
Index: trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScatteringModel.cc,v 1.39 2009/05/10 16:09:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4CoulombScatteringModel.cc,v 1.43 2009/10/28 10:14:13 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -39,4 +39,5 @@
 //
 // Modifications:
+//
 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
 //          logic of building - only elements from G4ElementTable
@@ -45,4 +46,6 @@
 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 
 // 09.06.08 V.Ivanchenko SelectIsotope is moved to the base class
+// 16.06.09 Consolandi rows 109, 111-112, 183, 185-186
+//
 //
 // Class Description:
@@ -84,6 +87,6 @@
 {
   SetupParticle(p);
-  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
-  SetupKinematic(ekin, cutEnergy);
+  if(kinEnergy < lowEnergyLimit) return 0.0;
+  SetupKinematic(kinEnergy, cutEnergy);
 
   // save lab system kinematics
@@ -103,5 +106,9 @@
   mom2 = momCM*momCM;
   tkin = sqrt(mom2 + m12) - mass;
-  invbeta2 = 1.0 +  m12/mom2;
+
+  //invbeta2 = 1.0 +  m12/mom2;
+ 
+  G4double fm = m2/(mass + m2);
+  invbeta2 = 1.0 +  m12*fm*fm/mom2;
 
   SetupTarget(Z, tkin);
@@ -127,12 +134,12 @@
 {
   G4double kinEnergy = dp->GetKineticEnergy();
-  if(kinEnergy <= DBL_MIN) return;
+  if(kinEnergy < lowEnergyLimit) return;
   DefineMaterial(couple);
   SetupParticle(dp->GetDefinition());
-  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
-  SetupKinematic(ekin, cutEnergy);
+  SetupKinematic(kinEnergy, cutEnergy);
 
   // Choose nucleus
-  currentElement = SelectRandomAtom(couple,particle,ekin,ecut,tkin);
+  currentElement = SelectRandomAtom(couple,particle,
+				    kinEnergy,ecut,kinEnergy);
 
   G4double Z  = currentElement->GetZ();
@@ -160,5 +167,4 @@
   G4double z1   = 1.0 - cost;
   if(z1 < 0.0) return;
-
   G4double sint = sqrt(z1*(1.0 + cost));
   G4double phi  = twopi * G4UniformRand();
@@ -175,15 +181,23 @@
   fParticleChange->ProposeMomentumDirection(newDirection);   
 
-  G4double elab = gam*(eCM + bet*pzCM);
-  ekin = elab - mass;
-  if(ekin < 0.0) ekin = 0.0;
-  fParticleChange->SetProposedKineticEnergy(ekin);
+  //  G4double elab = gam*(eCM + bet*pzCM);
+
+  G4double Ecm  = sqrt(mass*mass + m2*m2 + 2.0*etot*m2);
+  G4double elab = etot - m2*(ptot/Ecm)*(ptot/Ecm)*(1.-cost) ;
+
+ 
+  G4double finalT = elab - mass;
+  if(finalT < 0.0) finalT = 0.0;
+  fParticleChange->SetProposedKineticEnergy(finalT);
 
   // recoil
-  G4double erec = kinEnergy - ekin;
-
-  if(erec > recoilThreshold) {
+  G4double erec = kinEnergy - finalT;
+
+  G4double tcut = recoilThreshold;
+  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); } 
+
+  if(erec > tcut) {
     G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
-    G4double plab = sqrt(ekin*(ekin + 2.0*mass));
+    G4double plab = sqrt(finalT*(finalT + 2.0*mass));
     G4ThreeVector p2 = (ptot*dir - plab*newDirection).unit();
     G4DynamicParticle* newdp  = new G4DynamicParticle(ion, p2, erec);
Index: trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4GammaConversion.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
Index: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4GoudsmitSaundersonMscModel.cc,v 1.7 2009/04/20 19:22:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4GoudsmitSaundersonMscModel.cc,v 1.19 2009/11/09 18:40:25 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -44,6 +44,14 @@
 //                     step into two sub-steps occur only for msc regime
 //
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
+//                     adding a theta min limit due to screening effect of the atomic nucleus
+// 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
+//                     within CalculateIntegrals method
+// 05.10.2009 O.Kadri: tuning small angle theta distributions
+//                     assuming the case of lambdan<1 as single scattering regime
+//                     tuning theta sampling for theta below the screening angle
+//
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 
 //REFERENCES:
 //Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
@@ -62,5 +70,4 @@
 #include "G4MaterialCutsCouple.hh"
 #include "G4DynamicParticle.hh"
-#include "G4DataInterpolation.hh" 
 #include "G4Electron.hh"
 #include "G4Positron.hh"
@@ -82,7 +89,7 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
-  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(GeV),isInitialized(false)
+  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV),isInitialized(false)
 { 
-  fr=0.02,rangeinit=0.,masslimite=0.6*MeV;
+  fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
   particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
   tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
@@ -117,6 +124,7 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
-                          G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
+G4double 
+G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
+                       G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
 {  
   //Build cross section table : Taken from Ref.7
@@ -126,5 +134,4 @@
   if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
 
-  //value0=Lambda0;value1=Lambda1
   G4double value0,value1;
   CalculateIntegrals(p,Z,kinEnergy,value0,value1);
@@ -136,6 +143,7 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
-				       G4double safety)
+void 
+G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
+					       G4double safety)
 {    
   G4double kineticEnergy = dynParticle->GetKineticEnergy();
@@ -144,21 +152,19 @@
   G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
   G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
-  G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2,nu_interm;
+  G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2;
 
   ///////////////////////////////////////////
   // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4
   G4double  eloss = theManager->GetEnergy(particle,tPathLength,currentCouple);
-  if(eloss>0.5*kineticEnergy)return;
+  if(eloss>0.5*kineticEnergy)eloss=kineticEnergy-eloss;//exchange possibility between target atomic e- and incident particle
   G4double ee       = kineticEnergy - 0.5*eloss;
   G4double ttau     = ee/electron_mass_c2;
   G4double ttau2    = ttau*ttau;
   G4double epsilonpp= eloss/ee;
-  G4double temp2  = 0.166666*(4+ttau*(6+ttau*(7+ttau*(4+ttau))))*(epsilonpp/(ttau+1)/(ttau+2))*(epsilonpp/(ttau+1)/(ttau+2));
   G4double cst1=epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
 
   kineticEnergy *= (1 - cst1);
-  tPathLength *= (1 - temp2);
   ///////////////////////////////////////////
-  // additivity rule for mixture xsection calculation
+  // additivity rule for mixture and compound xsection calculation
   const G4Material* mat = currentCouple->GetMaterial();
   G4int nelm = mat->GetNumberOfElements();
@@ -166,6 +172,5 @@
   const G4double* theFraction = mat->GetFractionVector();
   G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume();
-  G4double scrA,llambda0,llambda1;
-  scrA=0.0;llambda0 =0.;llambda1=0.;
+  G4double llambda0=0.,llambda1=0.;
   for(G4int i=0;i<nelm;i++)
     {
@@ -175,60 +180,72 @@
       llambda1 += (theFraction[i]/l1);
     } 
-  if(llambda0>DBL_MIN)llambda0 =1./llambda0;
-  if(llambda1>DBL_MIN)llambda1 =1./llambda1;
-  G4double g1=llambda0/llambda1;
+  if(llambda0>DBL_MIN) llambda0 =1./llambda0;
+  if(llambda1>DBL_MIN) llambda1 =1./llambda1;
+  G4double g1=0.0;
+  if(llambda1>DBL_MIN) g1 = llambda0/llambda1;
+
   G4double x1,x0;
-
   x0=g1/2.;
   do
     {  
-      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)
+      x1 = x0-(x0*((1.+x0)*log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)
       delta = std::abs( x1 - x0 );    
       x0 = x1;  // new approximation becomes the old approximation for the next iteration
     } while (delta > 1e-10);
-  scrA = x1;
+  G4double scrA = x1;
 
   G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
   G4double lambdan=0.;
+  G4bool mscatt=false,noscatt=false;
+
   if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0;
-  if((lambdan<=1.0e-12)||(lambdan>1.0e+5))return;
-  G4bool noscatt=false;
-  G4bool singlescatt=false;
-  G4bool mscatt=false;
+  if((lambdan<=1.0e-12))return;
 
   G4double epsilon1=G4UniformRand();
-  if(epsilon1<(exp(-lambdan)))noscatt=true;// no scattering     
-  else if(epsilon1<((1.+lambdan)*exp(-lambdan)))//single scattering
-    {singlescatt=true;
-     ws=G4UniformRand();
-     ws= 1.-2.*scrA*ws/(1.-ws + scrA);
-     G4double phi0=twopi*G4UniformRand(); 
-     us=sqrt(1.-ws*ws)*cos(phi0);
-     vs=sqrt(1.-ws*ws)*sin(phi0);
-     G4double rr=G4UniformRand();
-     x_coord=(rr*us);
-     y_coord=(rr*vs);
-     z_coord=((1.-rr)+rr*ws);     
+  G4double expn = exp(-lambdan);
+  if((epsilon1<expn)||insideskin)// no scattering 
+    {noscatt=true;}
+  else if((epsilon1<((1.+lambdan)*expn)||(lambdan<1.)))
+    {
+      mscatt=false;
+      ws=G4UniformRand();
+      ws= 1.-2.*scrA*ws/(1.-ws + scrA);
+      if(acos(ws)<sqrt(scrA))//small angle approximation for theta less than screening angle
+      {G4int i=0;
+      do{i++;
+      ws=1.+0.5*atomPerVolume*tPathLength*log(G4UniformRand())/llambda1;
+      }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
+      if(i==19)ws=cos(sqrt(scrA));
+      }
+      G4double phi0=twopi*G4UniformRand(); 
+      us=sqrt(1.-ws*ws)*cos(phi0);
+      vs=sqrt(1.-ws*ws)*sin(phi0);
+      G4double rr=G4UniformRand();
+      x_coord=(rr*us);
+      y_coord=(rr*vs);
+      z_coord=((1.-rr)+rr*ws);
     }
   else
-  {mscatt=true;
-  // Ref.2 subsection 4.4 "The best solution found"
-  // Sample first substep scattering angle
-  SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
-  phi1  = twopi*G4UniformRand();
-  cosPhi1 = cos(phi1);
-  sinPhi1 = sin(phi1);
-
-  // Sample second substep scattering angle
-  SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
-  phi2  = twopi*G4UniformRand();
-  cosPhi2 = cos(phi2);
-  sinPhi2 = sin(phi2);
-
-  // Scattering direction
-  us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
-  vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
-  ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
-  }
+    {
+      mscatt=true;
+      // Ref.2 subsection 4.4 "The best solution found"
+      // Sample first substep scattering angle
+      SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
+      phi1  = twopi*G4UniformRand();
+      cosPhi1 = cos(phi1);
+      sinPhi1 = sin(phi1);
+
+      // Sample second substep scattering angle
+      SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
+      phi2  = twopi*G4UniformRand();
+      cosPhi2 = cos(phi2);
+      sinPhi2 = sin(phi2);
+
+      // Scattering direction
+      us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
+      vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
+      ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
+    }
+    
   G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
   G4ThreeVector newDirection(us,vs,ws);
@@ -238,37 +255,41 @@
   if((safety > tlimitminfix)&&(latDisplasment))
     {  
-      // Scattering coordinates
+
       if(mscatt)
-      {
-      if(scrA<DBL_MIN)scrA=DBL_MIN;
-      q1       = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
-      if(q1<DBL_MIN)q1=DBL_MIN;
-      Gamma    = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
-      Eta      = atomPerVolume*tPathLength/llambda1;      
-      delta    = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
-
-      nu = G4UniformRand(); 
-      nu = std::sqrt(nu);
-      nu0 = (1.0 - nu)/2.;
-      nu1 = nu*delta;
-      nu2 = nu*(1.0-delta);
-      nu_interm = 1.0 - nu0 - nu1 - nu2; 
-      x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu_interm*us);
-      y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu_interm*vs);
-      z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu_interm*ws)  ;
-      }
-      G4double r=sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
-
-      G4double check= 1.- tPathLength/zPathLength;
-      if(check<=0.)    return;
-      else if(r>check) {r=check;x_coord /=r;y_coord /=r;z_coord /=r;}
-
-      r *=tPathLength;
-      
+	{
+	  if(scrA<DBL_MIN)scrA=DBL_MIN;
+	  if(llambda1<DBL_MIN)llambda1=DBL_MIN;
+	  q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
+	  if(q1<DBL_MIN)q1=DBL_MIN;
+	  Gamma  = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
+	  Eta    = atomPerVolume*tPathLength/llambda1;      
+	  delta  = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
+
+	  nu = G4UniformRand(); 
+	  nu = std::sqrt(nu);
+	  nu0 = (1.0 - nu)/2.;
+	  nu1 = nu*delta;
+	  nu2 = nu*(1.0-delta);
+	  x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu0*us);
+	  y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu0*vs);
+	  z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu0*ws)  ;
+	}
+
+      // displacement is computed relatively to the end point
+      if(!noscatt)z_coord -= 1.0;//for noscatt zcoord z_coord !=0.
+      G4double rr = sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
+      G4double r  = rr*zPathLength;
+      /*
+      G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
+	     << " sinTheta= " << sqrt(1.0 - ws*ws) << " r(mm)= " << r
+	     << " trueStep(mm)= " << tPathLength
+	     << " geomStep(mm)= " << zPathLength
+	     << G4endl;
+      */
       if(r > tlimitminfix) {
 
-	G4ThreeVector latDirection = G4ThreeVector(x_coord*tPathLength,y_coord*tPathLength,z_coord*tPathLength);
-	latDirection.rotateUz(oldDirection);
- 
+        G4ThreeVector latDirection(x_coord/rr,y_coord/rr,z_coord/rr);
+        latDirection.rotateUz(oldDirection);
+
 	ComputeDisplacement(fParticleChange, latDirection, r, safety);
       }     
@@ -277,6 +298,8 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-void G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
-						     G4double &cost, G4double &sint)
+
+void 
+G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
+						G4double &cost, G4double &sint)
 {
   G4double u,Qn1,r1,tet;
@@ -284,26 +307,16 @@
   Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
 
-  if((lambdan<1.)||(Qn1<0.001))//plural scatt. or small angle scatt.
-	    {
-	      G4double xi1,lambdai=0.;
-	      G4int i=0;
-	      do {xi1=G4UniformRand();
-	      lambdai -=std::log(xi1);
-	      xi +=2.*scrA*xi1/(1.-xi1 + scrA);
-	      i++;
-	      }while((lambdai<lambdan)&&(i<30));
-	    
-	    }
-	  else {
-	    if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
-	    else{// procedure described by Benedito in Ref.1
-	         do{r1=G4UniformRand();
-	         u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
- 	         xi = 2.*u;
-       	         tet=acos(1.-xi);
-                 }while(tet*r1*r1>sin(tet));
-	       } 
-	      }
-
+if (Qn1<0.001)xi=-0.5*Qn1*log(G4UniformRand());
+else if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
+else
+{
+      // procedure described by Benedito in Ref.1
+      do{
+	r1=G4UniformRand();
+	u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
+	xi = 2.*u;
+	tet=acos(1.-xi);
+      }while(tet*r1*r1>sin(tet));
+}
 
   if(xi<0.)xi=0.;
@@ -315,53 +328,97 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-// Cubic spline log-log interpolation of Lambda0 and Lambda1
-// Screening parameter calculated according to Eq. 37 of Ref.1 
-void G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 
-						      G4double kinEnergy,G4double &Lam0,
-						      G4double &Lam1)
+// Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
+// linear log-log extrapolation between 1 GeV - 100 TeV
+
+void 
+G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 
+						 G4double kinEnergy,G4double &Lam0,
+						 G4double &Lam1)
 { 
-  //////// BEGIN OF: LAMBDA CALCULATION ////////////////////////////////
-  G4double TCSEForThisAtom[106],FTCSEForThisAtom[106],TCSPForThisAtom[106],FTCSPForThisAtom[106];
   G4double summ00=0.0;
   G4double summ10=0.0;
-  G4double InterpolatedValue=0.0;
-  
-  
-  G4int  iZ = G4int(Z);
-  if(iZ > 103) iZ = 103;
-  for(G4int i=0;i<106;i++)
-    {
-      TCSEForThisAtom[i]=TCSE[iZ-1][i];FTCSEForThisAtom[i]=FTCSE[iZ-1][i];
-      TCSPForThisAtom[i]=TCSP[iZ-1][i];FTCSPForThisAtom[i]=FTCSP[iZ-1][i];
-    }
-
+  G4double x1,x2,y1,y2,acoeff,bcoeff;
   G4double kineticE = kinEnergy;
   if(kineticE<lowKEnergy)kineticE=lowKEnergy;
   if(kineticE>highKEnergy)kineticE=highKEnergy;
   kineticE /= eV;
+  G4double logE=log(kineticE);
+  
+  G4int  iZ = G4int(Z);
+  if(iZ > 103) iZ = 103;
+
+  G4int enerInd=0;
+  for(G4int i=1;i<106;i++)
+  {
+  if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
+  }
+
+  if(p==G4Electron::Electron())        
+    {
+    if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b
+      {
+	x1=ener[enerInd];
+	x2=ener[enerInd+1];       
+	y1=TCSE[iZ-1][enerInd];
+	y2=TCSE[iZ-1][enerInd+1];
+	acoeff=(y2-y1)/(x2*x2-x1*x1);
+	bcoeff=y2-acoeff*x2*x2;
+	summ00=acoeff*logE*logE+bcoeff;
+	summ00 =exp(summ00);
+	y1=FTCSE[iZ-1][enerInd];
+	y2=FTCSE[iZ-1][enerInd+1];
+	acoeff=(y2-y1)/(x2*x2-x1*x1);
+	bcoeff=y2-acoeff*x2*x2;
+	summ10=acoeff*logE*logE+bcoeff;
+	summ10 =exp(summ10);
+      }
+    else  //Interpolation of the form y=ax+b
+      {  
+	x1=ener[104];
+	x2=ener[105];       
+	y1=TCSE[iZ-1][104];
+	y2=TCSE[iZ-1][105];
+	summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
+	summ00 =exp(summ00);
+	y1=FTCSE[iZ-1][104];
+	y2=FTCSE[iZ-1][105];
+	summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
+	summ10 =exp(summ10);
+      }
+    }
+  if(p==G4Positron::Positron())        
+    {
+    if(kinEnergy<=1.0e+9)
+      {
+	x1=ener[enerInd];
+	x2=ener[enerInd+1];       
+	y1=TCSP[iZ-1][enerInd];
+	y2=TCSP[iZ-1][enerInd+1];
+	acoeff=(y2-y1)/(x2*x2-x1*x1);
+	bcoeff=y2-acoeff*x2*x2;
+	summ00=acoeff*logE*logE+bcoeff;
+	summ00 =exp(summ00);
+	y1=FTCSP[iZ-1][enerInd];
+	y2=FTCSP[iZ-1][enerInd+1];
+	acoeff=(y2-y1)/(x2*x2-x1*x1);
+	bcoeff=y2-acoeff*x2*x2;
+	summ10=acoeff*logE*logE+bcoeff;
+	summ10 =exp(summ10);
+      }
+    else
+      {  
+	x1=ener[104];
+	x2=ener[105];       
+	y1=TCSP[iZ-1][104];
+	y2=TCSP[iZ-1][105];
+	summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
+	summ00 =exp(summ00);
+	y1=FTCSP[iZ-1][104];
+	y2=FTCSP[iZ-1][105];
+	summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
+	summ10 =exp(summ10);
+      }
+    }
     
-  if(p==G4Electron::Electron())        
-    {
-      MyValue= new G4DataInterpolation(ener,TCSEForThisAtom,106,0.0,0);  
-      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
-      delete  MyValue;
-      summ00 = std::exp(InterpolatedValue);  
-      MyValue= new G4DataInterpolation(ener,FTCSEForThisAtom,106,0.0,0);  
-      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
-      delete  MyValue;
-      summ10 = std::exp(InterpolatedValue);
-    }
-  if(p==G4Positron::Positron())        
-    {
-      MyValue= new G4DataInterpolation(ener,TCSPForThisAtom,106,0.0,0);  
-      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
-      delete  MyValue;
-      summ00 = std::exp(InterpolatedValue);  
-      MyValue= new G4DataInterpolation(ener,FTCSPForThisAtom,106,0.0,0);  
-      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
-      delete  MyValue;
-      summ10 = std::exp(InterpolatedValue);
-    }
-
   summ00 *=barn;
   summ10 *=barn;
@@ -374,7 +431,9 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //t->g->t step transformations taken from Ref.6 
-G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
-								  G4PhysicsTable* theTable,
-								  G4double currentMinimalStep)
+
+G4double 
+G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
+							 G4PhysicsTable* theTable,
+							 G4double currentMinimalStep)
 {
   tPathLength = currentMinimalStep;
@@ -406,4 +465,9 @@
 
   G4double presafety = sp->GetSafety();
+
+  //G4cout << "G4GS::StepLimit tPathLength= " 
+  //	 <<tPathLength<<" safety= " << presafety
+  //       << " range= " <<currentRange<< " lambda= "<<lambda1
+  //	 << " Alg: " << steppingAlgorithm <<G4endl;
 
   // far from geometry boundary
@@ -437,4 +501,17 @@
           else  smallstep = 1.;
 
+          //define stepmin here (it depends on lambda!)
+          //rough estimation of lambda_elastic/lambda_transport
+          G4double rat = currentKinEnergy/MeV ;
+          rat = 1.e-3/(rat*(10.+rat)) ;
+          //stepmin ~ lambda_elastic
+          stepmin = rat*lambda1;
+          skindepth = skin*stepmin;
+          //define tlimitmin
+          tlimitmin = 10.*stepmin;
+          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+
+	  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
+	  //	 << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
 	  // constraint from the geometry 
 	  if((geomlimit < geombig) && (geomlimit > geommin))
@@ -448,16 +525,4 @@
               tgeom = geombig;
 
-          //define stepmin here (it depends on lambda!)
-          //rough estimation of lambda_elastic/lambda_transport
-          G4double rat = currentKinEnergy/MeV ;
-          rat = 1.e-3/(rat*(10.+rat)) ;
-          //stepmin ~ lambda_elastic
-          stepmin = rat*lambda1;
-          skindepth = skin*stepmin;
-
-          //define tlimitmin
-          tlimitmin = 10.*stepmin;
-          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
-
         }
 
@@ -471,4 +536,7 @@
 
       if(tlimit > tgeom) tlimit = tgeom;
+
+      //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
+      //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
 
       // shortcut
@@ -499,5 +567,5 @@
       if(tlimit < stepmin) tlimit = stepmin;
 
-      if(tPathLength > tlimit) tPathLength = tlimit  ; 
+      if(tPathLength > tlimit) tPathLength = tlimit; 
 
     }
@@ -519,22 +587,22 @@
 
       if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
-      {
-        rangeinit = currentRange;
-        fr = facrange;
-        // 9.1 like stepping for e+/e- only (not for muons,hadrons)
-        if(mass < masslimite) 
-        {
-          if(lambda1 > currentRange)
-            rangeinit = lambda1;
-          if(lambda1 > lambdalimit)
-            fr *= 0.75+0.25*lambda1/lambdalimit;
-        }
-
-        //lower limit for tlimit
-        G4double rat = currentKinEnergy/MeV ;
-        rat = 1.e-3/(rat*(10.+rat)) ;
-        tlimitmin = 10.*lambda1*rat;
-        if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
-      }
+	{
+	  rangeinit = currentRange;
+	  fr = facrange;
+	  // 9.1 like stepping for e+/e- only (not for muons,hadrons)
+	  if(mass < masslimite) 
+	    {
+	      if(lambda1 > currentRange)
+		rangeinit = lambda1;
+	      if(lambda1 > lambdalimit)
+		fr *= 0.75+0.25*lambda1/lambdalimit;
+	    }
+
+	  //lower limit for tlimit
+	  G4double rat = currentKinEnergy/MeV ;
+	  rat = 1.e-3/(rat*(10.+rat)) ;
+	  tlimitmin = 10.*lambda1*rat;
+	  if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+	}
       //step limit
       tlimit = fr*rangeinit;               
@@ -561,8 +629,11 @@
 	}
     }
+  //G4cout << "tPathLength= " << tPathLength  
+  // << " currentMinimalStep= " << currentMinimalStep << G4endl;
   return tPathLength ;
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
 {
@@ -602,5 +673,6 @@
       zmean = 1./(par1*par3) ;
   } else {
-    G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,currentCouple);
+    G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,
+					currentCouple);
 
     lambda11 = GetLambda(T1);
@@ -614,14 +686,14 @@
   zPathLength = zmean ;
   //  sample z
-  if(samplez)
-  {
+  if(samplez) {
+
     const G4double  ztmax = 0.99, onethird = 1./3. ;
     G4double zt = zmean/tPathLength ;
 
-    if (tPathLength > stepmin && zt < ztmax)              
-    {
+    if (tPathLength > stepmin && zt < ztmax) {
+
       G4double u,cz1;
-      if(zt >= onethird)
-      {
+      if(zt >= onethird) {
+
         G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
         cz1 = 1.+cz ;
@@ -629,10 +701,9 @@
         G4double grej ;
         do {
-            u = exp(log(G4UniformRand())/cz1) ;
-            grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
-           } while (grej < G4UniformRand()) ;
-      }
-      else
-      {
+	  u = exp(log(G4UniformRand())/cz1) ;
+	  grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
+	} while (grej < G4UniformRand()) ;
+
+      } else {
         cz1 = 1./zt-1.;
         u = 1.-exp(log(G4UniformRand())/cz1) ;
@@ -642,11 +713,13 @@
   }
   if(zPathLength > lambda1) zPathLength = lambda1;
-
- return zPathLength;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
+  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
+
+  return zPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double 
+G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
 {
   // step defined other than transportation 
@@ -673,4 +746,5 @@
     }
   if(tPathLength < geomStepLength) tPathLength = geomStepLength;
+  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
 
   return tPathLength;
@@ -678,72 +752,71 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//Total & first transport x sections for e-/e+ generated from ELSEPA code
+
 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
 { 
-  ///////////////////////////////////////
-  //Total & first transport x sections of e-/e+ from ELSEPA code
   G4String filename = "XSECTIONS.dat";
 
-    char* path = getenv("G4LEDATA");
-    if (!path)
-      {
-        G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set";
-        G4Exception(excep);
-      }
-
-    G4String pathString(path);
-    G4String dirFile = pathString + "/msc_GS/" + filename;
-    FILE *infile;
-    infile = fopen(dirFile,"r"); 
-    if (infile == 0)
-      {
-	G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
-	G4Exception(excep);
-      }
-
-    // Read parameters from tables and take logarithms
-    G4float aRead;
-    for(G4int i=0 ; i<106 ;i++){
-	  fscanf(infile,"%f\t",&aRead);
-          if(aRead > 0.0) aRead = std::log(aRead);
-          else  aRead = 0.0;
-          ener[i]=aRead;
+  char* path = getenv("G4LEDATA");
+  if (!path)
+    {
+      G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set properly";
+      G4Exception(excep);
+    }
+
+  G4String pathString(path);
+  G4String dirFile = pathString + "/msc_GS/" + filename;
+  FILE *infile;
+  infile = fopen(dirFile,"r"); 
+  if (infile == 0)
+    {
+      G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
+      G4Exception(excep);
+    }
+
+  // Read parameters from tables and take logarithms
+  G4float aRead;
+  for(G4int i=0 ; i<106 ;i++){
+    fscanf(infile,"%f\t",&aRead);
+    if(aRead > 0.0) aRead = log(aRead);
+    else  aRead = 0.0;
+    ener[i]=aRead;
+  }        
+  for(G4int j=0;j<103;j++){
+    for(G4int i=0;i<106;i++){
+      fscanf(infile,"%f\t",&aRead);
+      if(aRead > 0.0) aRead = log(aRead);
+      else  aRead = 0.0;
+      TCSE[j][i]=aRead;
     }        
-    for(G4int j=0;j<103;j++){
-      for(G4int i=0;i<106;i++){
-	  fscanf(infile,"%f\t",&aRead);
-          if(aRead > 0.0) aRead = std::log(aRead);
-          else  aRead = 0.0;
-	  TCSE[j][i]=aRead;
-	}        
-     }
-    for(G4int j=0;j<103;j++){
-      for(G4int i=0;i<106;i++){
-	  fscanf(infile,"%f\t",&aRead);
-          if(aRead > 0.0) aRead = std::log(aRead);
-          else  aRead = 0.0;
-	  FTCSE[j][i]=aRead;      
-	}        
-      }    
-    for(G4int j=0;j<103;j++){
-      for(G4int i=0;i<106;i++){
-	  fscanf(infile,"%f\t",&aRead);
-          if(aRead > 0.0) aRead = std::log(aRead);
-          else  aRead = 0.0;
-	  TCSP[j][i]=aRead;      
-	}        
-     }
-    for(G4int j=0;j<103;j++){
-      for(G4int i=0;i<106;i++){
-	  fscanf(infile,"%f\t",&aRead);
-          if(aRead > 0.0) aRead = std::log(aRead);
-          else  aRead = 0.0;
-	  FTCSP[j][i]=aRead;      
-	}        
-     }
-
-    fclose(infile);
-   //End loading XSections and Energies
-
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+  }
+  for(G4int j=0;j<103;j++){
+    for(G4int i=0;i<106;i++){
+      fscanf(infile,"%f\t",&aRead);
+      if(aRead > 0.0) aRead = log(aRead);
+      else  aRead = 0.0;
+      FTCSE[j][i]=aRead;      
+    }        
+  }    
+  for(G4int j=0;j<103;j++){
+    for(G4int i=0;i<106;i++){
+      fscanf(infile,"%f\t",&aRead);
+      if(aRead > 0.0) aRead = log(aRead);
+      else  aRead = 0.0;
+      TCSP[j][i]=aRead;      
+    }        
+  }
+  for(G4int j=0;j<103;j++){
+    for(G4int i=0;i<106;i++){
+      fscanf(infile,"%f\t",&aRead);
+      if(aRead > 0.0) aRead = log(aRead);
+      else  aRead = 0.0;
+      FTCSP[j][i]=aRead;      
+    }        
+  }
+
+  fclose(infile);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4GoudsmitSaundersonTable.cc,v 1.1 2009/03/05 18:48:30 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4GoudsmitSaundersonTable.cc,v 1.4 2009/08/28 16:36:52 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -39,4 +39,6 @@
 // Modifications:
 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
+// 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root finding parameter error's 
+//                     within SampleTheta method
 //
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -167,12 +169,16 @@
   0.996273,0.997261,0.998097,0.998782,0.999315,0.999695,0.999924,1.0};
 
-G4double G4GoudsmitSaundersonTable::PDF[]  = {-1.};
-G4double G4GoudsmitSaundersonTable::CPDF[] = {-1.};
+G4double* G4GoudsmitSaundersonTable::PDF  = 0;
+G4double* G4GoudsmitSaundersonTable::CPDF = 0;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable()
 {
-  if(PDF[0] < 0.0){ 
+  if(PDF == 0){ 
     G4cout << "### G4GoudsmitSaundersonTable loading PDF data" << G4endl;
+
+    PDF = new G4double [76*11*320];
+    CPDF= new G4double [76*11*320];
+
     LoadPDFandCPDFdata();
   }
@@ -181,5 +187,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable()
-{}
+{
+  if(PDF) {
+    delete [] PDF;
+    delete [] CPDF;
+    PDF = 0;
+  }
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -194,7 +206,7 @@
   ///////////////////////////////////////////////////////////////////////////
   // Find Lambda and Chia2 Index
-  for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1]))KIndex=k;}  
+  for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1])){KIndex=k;break;}}  
   for(G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];}  
-  for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1]))JIndex=j;}
+  for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1])){JIndex=j;break;}}
 
   ///////////////////////////////////////////////////////////////////////////
@@ -218,12 +230,13 @@
     ThisPDF[i]=Ckj*PDF1[i]+CkjPlus1*PDF2[i]+CkPlus1j*PDF3[i]+CkPlus1jPlus1*PDF4[i];
     ThisCPDF[i]=Ckj*CPDF1[i]+CkjPlus1*CPDF2[i]+CkPlus1j*CPDF3[i]+CkPlus1jPlus1*CPDF4[i];  
+  // Find u Index using secant method
+  if((i!=0)&&((rndm>=ThisCPDF[i-1])&&(rndm<ThisCPDF[i])))  {IIndex=i-1;break;}
   }
 
   ///////////////////////////////////////////////////////////////////////////
-  // Find u Index using secant method
-  for(G4int i=0;i<319;i++){if((rndm>=ThisCPDF[i])&&(rndm<ThisCPDF[i+1]))IIndex=i;}
   //CPDF^-1(rndm)=x ==> CPDF(x)=rndm;
   aa=uvalues[IIndex];
   b=uvalues[IIndex+1];
+  
   do{
     m=0.5*(aa+b);
@@ -233,5 +246,5 @@
     if(F>0.)b=m;
     else aa=m;
-  } while(sqrt((b-aa)*(b-aa))>1.0e-10);
+  } while(sqrt((b-aa)*(b-aa))>1.0e-6);
 
   return m;
Index: trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4HeatedKleinNishinaCompton.cc,v 1.5 2009/04/12 17:09:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4InitXscPAI.cc	(revision 1196)
@@ -26,5 +26,5 @@
 //
 // $Id: G4InitXscPAI.cc,v 1.9 2006/06/29 19:53:00 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
Index: trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4IonFluctuations.cc,v 1.26 2009/03/31 13:24:40 toshito Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4KleinNishinaCompton.cc,v 1.10 2009/05/15 17:12:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MollerBhabhaModel.cc,v 1.34 2009/04/24 17:15:46 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4MollerBhabhaModel.cc,v 1.35 2009/11/09 19:16:13 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -258,15 +258,16 @@
 
   //density correction 
-  G4double cden  = material->GetIonisation()->GetCdensity();
-  G4double mden  = material->GetIonisation()->GetMdensity();
-  G4double aden  = material->GetIonisation()->GetAdensity();
-  G4double x0den = material->GetIonisation()->GetX0density();
-  G4double x1den = material->GetIonisation()->GetX1density();
+  //G4double cden  = material->GetIonisation()->GetCdensity();
+  //G4double mden  = material->GetIonisation()->GetMdensity();
+  //G4double aden  = material->GetIonisation()->GetAdensity();
+  //G4double x0den = material->GetIonisation()->GetX0density();
+  //G4double x1den = material->GetIonisation()->GetX1density();
   G4double x     = log(bg2)/twoln10;
 
-  if (x >= x0den) {
-    dedx -= twoln10*x - cden;
-    if (x < x1den) dedx -= aden*pow(x1den-x, mden);
-  } 
+  //if (x >= x0den) {
+  //  dedx -= twoln10*x - cden;
+  //  if (x < x1den) dedx -= aden*pow(x1den-x, mden);
+  //}
+  dedx -= material->GetIonisation()->DensityCorrection(x); 
 
   // now you can compute the total ionization loss
Index: trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MscModel71.cc,v 1.7 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4MscModel71.cc,v 1.8 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -104,5 +104,5 @@
 G4MscModel71::G4MscModel71(G4double& m_dtrl, G4double& m_NuclCorrPar,
 			   G4double& m_FactPar, G4double& m_factail,
-		       G4bool& m_samplez, const G4String& nam)
+			   G4bool& m_samplez, const G4String& nam)
   : G4VEmModel(nam),
   taubig(8.0),
@@ -117,5 +117,9 @@
 {
   stepmin       = 1.e-6*mm;
-  currentRange  = 0. ;
+  currentRange  = 0.;
+  G4cout << G4endl;
+  G4cout << "!!! G4MscModel71 class is obsolete and will be removed for the next major Geant4 release !!!" << G4endl;
+  G4cout << "!!! Please use other models (G4UrbanMscModel90, 92, 93) !!!" << G4endl;
+  G4cout << G4endl;
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MultipleScattering.cc,v 1.75 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4MultipleScattering.cc,v 1.77 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -----------------------------------------------------------------------------
@@ -129,5 +129,5 @@
 
 #include "G4MultipleScattering.hh"
-#include "G4UrbanMscModel.hh"
+#include "G4UrbanMscModel92.hh"
 #include "G4MscStepLimitType.hh"
 #include "G4UrbanMscModel.hh"
@@ -141,4 +141,10 @@
 {
   isInitialized = false;  
+  G4cout << G4endl;
+  G4cout << "!!! G4MultipleScattering class is obsolete and will be removed for the next major Geant4 release !!!" << G4endl;
+  G4cout << "!!! Please use G4eMultipleScattering for e+ and e- !!!" << G4endl;
+  G4cout << "!!! Please use G4MuMultipleScattering for mu+ and mu- !!!" << G4endl;
+  G4cout << "!!! Please use G4hMultipleScattering for hadrons and ions !!!" << G4endl;
+  G4cout << G4endl;
 }
 
@@ -157,34 +163,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4MultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)
+void G4MultipleScattering::InitialiseProcess(const G4ParticleDefinition*)
 {
-  // Modification of parameters between runs
-  if(isInitialized) {
-    if (p->GetParticleType() != "nucleus") {
-      mscUrban->SetStepLimitType(StepLimitType());
-      mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-      mscUrban->SetSkin(Skin());
-      mscUrban->SetRangeFactor(RangeFactor());
-      mscUrban->SetGeomFactor(GeomFactor());
-    }
-    return;
-  }
-
-  // defaults for ions, which cannot be overwritten
-  if (p->GetParticleType() == "nucleus") {
-    SetStepLimitType(fMinimal);
-    SetLateralDisplasmentFlag(false);
-    SetBuildLambdaTable(false);
-  }
+  if(isInitialized) return;
 
   // initialisation of parameters - defaults for particles other
   // than ions can be overwritten by users
-  mscUrban = new G4UrbanMscModel();
-  mscUrban->SetStepLimitType(StepLimitType());
-  mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-  mscUrban->SetSkin(Skin());
-  mscUrban->SetRangeFactor(RangeFactor());
-  mscUrban->SetGeomFactor(GeomFactor());
-
+  mscUrban = new G4UrbanMscModel92();
   AddEmModel(1,mscUrban);
   isInitialized = true;
Index: trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering71.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering71.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4MultipleScattering71.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4MultipleScattering71.cc,v 1.5 2008/07/16 11:27:41 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -----------------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4PAIModel.cc,v 1.47 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4PAIModel.cc,v 1.51 2009/08/12 21:28:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -41,4 +41,5 @@
 // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary
 // 08.04.05 Major optimisation of internal interfaces (V.Ivantchenko)
+// 26.07.09 Fixed logic to work with several materials (V.Ivantchenko)
 //
 
@@ -63,5 +64,4 @@
 #include "G4ParticleDefinition.hh"
 #include "G4ParticleChangeForLoss.hh"
-#include "G4GeometryTolerance.hh"
 
 ////////////////////////////////////////////////////////////////////////
@@ -162,8 +162,4 @@
 
   // Prepare initialization
-
-  fPAItransferTable = new G4PhysicsTable(fTotBin);
-  fPAIdEdxTable = new G4PhysicsTable(fTotBin);
-
   const G4ProductionCutsTable* theCoupleTable =
         G4ProductionCutsTable::GetProductionCutsTable();
@@ -182,7 +178,10 @@
       //G4cout << "Reg <" <<curReg->GetName() << ">  mat <" 
       //	     << fMaterial->GetName() << ">  fCouple= " 
-      //	     << fCutCouple<<G4endl;
+      //	     << fCutCouple<<"  " << p->GetParticleName() <<G4endl;
       if( fCutCouple ) {
 	fMaterialCutsCoupleVector.push_back(fCutCouple);
+
+	fPAItransferTable = new G4PhysicsTable(fTotBin+1);
+	fPAIdEdxTable = new G4PhysicsTable(fTotBin+1);
 
 	fDeltaCutInKinEnergy = 
@@ -263,5 +262,4 @@
   G4double tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ;
 
-  if(fdEdxVector) delete fdEdxVector;
   fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy,
 					fHighestKineticEnergy,
@@ -272,5 +270,5 @@
   deltaLow = 100.*eV; // 0.5*eV ;
 
-  for (G4int i = 0 ; i < fTotBin ; i++)  //The loop for the kinetic energy
+  for (G4int i = 0 ; i <= fTotBin ; i++)  //The loop for the kinetic energy
   {
     LowEdgeEnergy = fParticleEnergyVector->GetLowEdgeEnergy(i) ;
@@ -283,17 +281,7 @@
     Tkin = Tmax ;
 
-    // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV"
-    // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl;
-  
     if ( Tmax < Tmin + deltaLow )  // low energy safety
       Tkin = Tmin + deltaLow ;
 
-    /*
-    G4PAIxSection protonPAI( fMatIndex,
-                             Tkin,
-                             bg2,
-                             fSandiaPhotoAbsCof,
-                             fSandiaIntervalNumber  ) ;
-    */
     fPAIySection.Initialize(fMaterial, Tkin, bg2);
 
@@ -325,8 +313,4 @@
 
   }                                        // end of Tkin loop
-  //  theLossTable->insert(fdEdxVector);
-  // end of material loop
-  // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ;
-  // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ;
 }
 
@@ -339,10 +323,4 @@
 void G4PAIModel::BuildLambdaVector()
 {
-  //G4double kCarTolerance = G4GeometryTolerance::GetInstance()
-  //                         ->GetSurfaceTolerance();
-
-  if (fLambdaVector)   delete fLambdaVector;
-  if (fdNdxCutVector)  delete fdNdxCutVector;
-
   fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy,
 					  fHighestKineticEnergy,
@@ -356,10 +334,12 @@
 	  <<fDeltaCutInKinEnergy/keV<<" keV"<<G4endl;
   }
-  for (G4int i = 0 ; i < fTotBin ; i++ )
+  for (G4int i = 0 ; i <= fTotBin ; i++ )
   {
     G4double dNdxCut = GetdNdxCut(i,fDeltaCutInKinEnergy) ;
-    G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
-
-    //    if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance ; // Mmm ???
+    //G4cout << "i= " << i << " x= " << dNdxCut << G4endl;
+    if(dNdxCut < 0.0) dNdxCut = 0.0; 
+    //    G4double lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut ;
+    G4double lambda = DBL_MAX;
+    if(dNdxCut > 0.0) lambda = 1.0/dNdxCut;
 
     fLambdaVector->PutValue(i, lambda) ;
@@ -395,4 +375,6 @@
   y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ;
   y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ;
+  if(y1 < 0.0) y1 = 0.0;
+  if(y2 < 0.0) y2 = 0.0;
   // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
   x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
@@ -409,4 +391,5 @@
   }
   //  G4cout<<""<<dNdxCut<<G4endl;
+  if(dNdxCut < 0.0) dNdxCut = 0.0; 
   return dNdxCut ;
 }
@@ -420,5 +403,5 @@
   G4int iTransfer;
   G4double x1, x2, y1, y2, dEdxCut;
-  // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
+  //G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl;
   // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength())
   //           <<G4endl;  
@@ -438,8 +421,10 @@
   y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ;
   y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ;
-  // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
+  if(y1 < 0.0) y1 = 0.0;
+  if(y2 < 0.0) y2 = 0.0;
+  //G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl;
   x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
   x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
-  // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
+  //G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl;
 
   if ( y1 == y2 )    dEdxCut = y2 ;
@@ -451,5 +436,6 @@
     else             dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
-  //  G4cout<<""<<dEdxCut<<G4endl;
+  //G4cout<<""<<dEdxCut<<G4endl;
+  if(dEdxCut < 0.0) dEdxCut = 0.0; 
   return dEdxCut ;
 }
@@ -463,5 +449,4 @@
 {
   G4int iTkin,iPlace;
-  size_t jMat;
   
   //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
@@ -474,18 +459,27 @@
   const G4MaterialCutsCouple* matCC = CurrentCouple();
 
-  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
+  size_t jMat = 0;
+  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
   {
     if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
   }
-  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
+  //G4cout << "jMat= " << jMat 
+  //	 << " jMax= " << fMaterialCutsCoupleVector.size()
+  //       << " matCC: " << matCC;
+  //if(matCC) G4cout << " mat: " << matCC->GetMaterial()->GetName();
+  //  G4cout << G4endl;
+  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
 
   fPAIdEdxTable = fPAIdEdxBank[jMat];
   fdEdxVector = fdEdxTable[jMat];
-  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
+  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
   {
     if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;    
   }
+  //G4cout << "E= " << scaledTkin << " iTkin= " << iTkin 
+  //	 << "  " << fdEdxVector->GetVectorLength()<<G4endl; 
   iPlace = iTkin - 1;
   if(iPlace < 0) iPlace = 0;
+  else if(iPlace >= fTotBin) iPlace = fTotBin-1;
   G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
   if( dEdx < 0.) dEdx = 0.;
@@ -502,5 +496,4 @@
 {
   G4int iTkin,iPlace;
-  size_t jMat;
   G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
   if(tmax <= cutEnergy) return 0.0;
@@ -511,13 +504,14 @@
   const G4MaterialCutsCouple* matCC = CurrentCouple();
 
-  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
+  size_t jMat = 0;
+  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
   {
     if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
   }
-  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
+  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
 
   fPAItransferTable = fPAIxscBank[jMat];
 
-  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
+  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
   {
     if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin)) break ;    
@@ -525,14 +519,15 @@
   iPlace = iTkin - 1;
   if(iPlace < 0) iPlace = 0;
-
-  // G4cout<<"iPlace = "<<iPlace<<"; tmax = "
+  else if(iPlace >= fTotBin) iPlace = fTotBin-1; 
+
+  //G4cout<<"iPlace = "<<iPlace<<"; tmax = "
   // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl;  
   cross1 = GetdNdxCut(iPlace,tmax) ;
-  // G4cout<<"cross1 = "<<cross1<<G4endl;  
+  //G4cout<<"cross1 = "<<cross1<<G4endl;  
   cross2 = GetdNdxCut(iPlace,cutEnergy) ;  
-  // G4cout<<"cross2 = "<<cross2<<G4endl;  
+  //G4cout<<"cross2 = "<<cross2<<G4endl;  
   cross  = (cross2-cross1)*charge2;
-  // G4cout<<"cross = "<<cross<<G4endl;  
-  if( cross < DBL_MIN) cross = DBL_MIN;
+  //G4cout<<"cross = "<<cross<<G4endl;  
+  if( cross < DBL_MIN) cross = 0.0;
   //  if( cross2 < DBL_MIN) cross2 = DBL_MIN;
 
@@ -552,10 +547,10 @@
 				   G4double maxEnergy)
 {
-  size_t jMat;
-  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
+  size_t jMat = 0;
+  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
   {
     if( matCC == fMaterialCutsCoupleVector[jMat] ) break;
   }
-  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
+  if(jMat == fMaterialCutsCoupleVector.size()) return;
 
   fPAItransferTable = fPAIxscBank[jMat];
@@ -636,5 +631,5 @@
   G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ;
 
-  for(iTkin=0;iTkin<fTotBin;iTkin++)
+  for(iTkin=0;iTkin<=fTotBin;iTkin++)
   {
     if(scaledTkin < fParticleEnergyVector->GetLowEdgeEnergy(iTkin))  break ;
@@ -643,4 +638,5 @@
   // G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
   if(iPlace < 0) iPlace = 0;
+  else if(iPlace >= fTotBin) iPlace = fTotBin-1; 
   dNdxCut1 = (*fdNdxCutVector)(iPlace) ;  
   // G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
@@ -758,20 +754,19 @@
                                                G4double&)
 {
-  size_t jMat;
-  for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
+  size_t jMat = 0;
+  for(;jMat < fMaterialCutsCoupleVector.size(); ++jMat )
   {
     if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break;
   }
-  if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--;
+  if(jMat == fMaterialCutsCoupleVector.size()) return 0.0;
 
   fPAItransferTable = fPAIxscBank[jMat];
   fdNdxCutVector   = fdNdxCutTable[jMat];
 
-  G4int iTkin, iTransfer, iPlace  ;
+  G4int iTkin, iTransfer, iPlace;
   G4long numOfCollisions=0;
 
-  //  G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
-  //G4cout<<"in:  "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl ;
-
+  //G4cout<<"G4PAIModel::SampleFluctuations"<<G4endl ;
+  //G4cout<<"in: "<<fMaterialCutsCoupleVector[jMat]->GetMaterial()->GetName()<<G4endl;
   G4double loss = 0.0, charge2 ;
   G4double stepSum = 0., stepDelta, lambda, omega; 
@@ -784,5 +779,5 @@
   G4double TkinScaled = Tkin*MassRatio ;
 
-  for(iTkin=0;iTkin<fTotBin;iTkin++)
+  for(iTkin=0;iTkin<=fTotBin;iTkin++)
   {
     if(TkinScaled < fParticleEnergyVector->GetLowEdgeEnergy(iTkin))   break ;
@@ -790,8 +785,8 @@
   iPlace = iTkin - 1 ; 
   if(iPlace < 0) iPlace = 0;
-  //  G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
+  else if(iPlace >= fTotBin) iPlace = fTotBin - 1;
+  //G4cout<<"from search, iPlace = "<<iPlace<<G4endl ;
   dNdxCut1 = (*fdNdxCutVector)(iPlace) ;  
-  //  G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
-
+  //G4cout<<"dNdxCut1 = "<<dNdxCut1<<G4endl ;
 
   if(iTkin == fTotBin) // Fermi plato, try from left
@@ -810,5 +805,5 @@
      numOfCollisions++;
     }   
-    //    G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
+    //G4cout<<"##1 numOfCollisions = "<<numOfCollisions<<G4endl ;
 
     while(numOfCollisions)
@@ -823,5 +818,5 @@
       }
       omega = GetEnergyTransfer(iPlace,position,iTransfer);
-      // G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
+      //G4cout<<"G4PAIModel::SampleFluctuations, omega = "<<omega/keV<<" keV; "<<"\t";
       loss += omega;
       numOfCollisions-- ;
@@ -831,5 +826,5 @@
   {
     dNdxCut2 = (*fdNdxCutVector)(iPlace+1) ; 
-    //  G4cout<<"dNdxCut2 = "<<dNdxCut2<<G4endl ;
+    //G4cout<<"dNdxCut2 = "<<dNdxCut2<< " iTkin= "<<iTkin<<" iPlace= "<<iPlace<<G4endl;
  
     if(iTkin == 0) // Tkin is too small, trying from right only
@@ -839,15 +834,15 @@
       //  numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ;
       //  numOfCollisions = G4Poisson(meanNumber) ;
-    if( meanNumber > 0.) lambda = step/meanNumber;
-    else                 lambda = DBL_MAX;
-    while(numb)
-    {
-     stepDelta = CLHEP::RandExponential::shoot(lambda);
-     stepSum += stepDelta;
-     if(stepSum >= step) break;
-     numOfCollisions++;
-    }   
-
-    //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
+      if( meanNumber > 0.) lambda = step/meanNumber;
+      else                 lambda = DBL_MAX;
+      while(numb)
+	{
+	  stepDelta = CLHEP::RandExponential::shoot(lambda);
+	  stepSum += stepDelta;
+	  if(stepSum >= step) break;
+	  numOfCollisions++;
+	}   
+
+      //G4cout<<"##2 numOfCollisions = "<<numOfCollisions<<G4endl ;
 
       while(numOfCollisions)
@@ -862,5 +857,5 @@
         }
         omega = GetEnergyTransfer(iPlace,position,iTransfer);
-        // G4cout<<omega/keV<<"\t";
+        //G4cout<<omega/keV<<"\t";
         loss += omega;
         numOfCollisions-- ;
@@ -875,7 +870,8 @@
       W2 = (TkinScaled - E1)*W ;
 
-      // G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
+      //G4cout << fPAItransferTable->size() << G4endl;
+      //G4cout<<"(*(*fPAItransferTable)(iPlace))(0) = "<<
       //   (*(*fPAItransferTable)(iPlace))(0)<<G4endl ;
-      // G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
+      //G4cout<<"(*(*fPAItransferTable)(iPlace+1))(0) = "<<
       //     (*(*fPAItransferTable)(iPlace+1))(0)<<G4endl ;
 
@@ -885,15 +881,15 @@
       //  numOfCollisions = RandPoisson::shoot(meanNumber) ;
       // numOfCollisions = G4Poisson(meanNumber) ;
-    if( meanNumber > 0.) lambda = step/meanNumber;
-    else                 lambda = DBL_MAX;
-    while(numb)
-    {
-     stepDelta = CLHEP::RandExponential::shoot(lambda);
-     stepSum += stepDelta;
-     if(stepSum >= step) break;
-     numOfCollisions++;
-    }   
-
-    //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
+      if( meanNumber > 0.) lambda = step/meanNumber;
+      else                 lambda = DBL_MAX;
+      while(numb)
+	{
+	  stepDelta = CLHEP::RandExponential::shoot(lambda);
+	  stepSum += stepDelta;
+	  if(stepSum >= step) break;
+	  numOfCollisions++;
+	}   
+
+      //G4cout<<"##3 numOfCollisions = "<<numOfCollisions<<endl ;
 
       while(numOfCollisions)
@@ -923,5 +919,5 @@
     }
   } 
-  //  G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
+  //G4cout<<"PAIModel AlongStepLoss = "<<loss/keV<<" keV, on step = "
   //  <<step/mm<<" mm"<<G4endl ; 
   if(loss > Tkin) loss=Tkin;
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4PAIPhotonModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4PAIPhotonModel.cc,v 1.23 2009/07/26 15:51:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -310,5 +310,5 @@
   deltaLow = 100.*eV; // 0.5*eV ;
 
-  for (G4int i = 0 ; i < fTotBin ; i++)  //The loop for the kinetic energy
+  for (G4int i = 0 ; i <= fTotBin ; i++)  //The loop for the kinetic energy
   {
     LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ;
@@ -441,5 +441,5 @@
 	  <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl;
   }
-  for ( i = 0 ; i < fTotBin ; i++ )
+  for ( i = 0 ; i <= fTotBin ; i++ )
   {
     dNdxPhotonCut  = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow);
@@ -664,5 +664,5 @@
   fPAIdEdxTable = fPAIdEdxBank[jMat];
   fdEdxVector = fdEdxTable[jMat];
-  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
+  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
   {
     if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;    
@@ -722,5 +722,5 @@
   fPAIplasmonTable  = fPAIplasmonBank[jMat];
 
-  for(iTkin = 0 ; iTkin < fTotBin ; iTkin++)
+  for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++)
   {
     if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ;    
@@ -789,5 +789,5 @@
 
   G4int iTkin;
-  for(iTkin=0;iTkin<fTotBin;iTkin++)
+  for(iTkin=0;iTkin<=fTotBin;iTkin++)
   {
     if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))  break ;
@@ -1052,5 +1052,5 @@
   G4double cof        = step*charge2;
 
-  for( iTkin = 0; iTkin < fTotBin; iTkin++)
+  for( iTkin = 0; iTkin <= fTotBin; iTkin++)
   {
     if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))   break ;
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIxSection.cc	(revision 1196)
@@ -26,5 +26,5 @@
 //
 // $Id: G4PAIxSection.cc,v 1.24 2008/05/30 16:04:40 grichine Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIySection.cc	(revision 1196)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4PAIySection.cc,v 1.3 2007/10/01 18:38:10 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4PAIySection.cc,v 1.4 2009/07/26 15:51:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // 
@@ -41,4 +41,6 @@
 //
 // 01.10.07, V.Ivanchenko create using V.Grichine G4PAIxSection class
+// 26.07.09, V.Ivanchenko added protection for mumerical exceptions for 
+//              low-density materials
 //
 
@@ -128,10 +130,6 @@
 	// G4cout<<fEnergyInterval[i]<<"\t"<<fA1[i]<<"\t"<<fA2[i]<<"\t"
 	//   <<fA3[i]<<"\t"<<fA4[i]<<"\t"<<G4endl ;
-        if(fEnergyInterval[i+1]-fEnergyInterval[i] >
+        if(fEnergyInterval[i+1]-fEnergyInterval[i] <
            1.5*fDelta*(fEnergyInterval[i+1]+fEnergyInterval[i]))
-	  {
-	    continue ;
-	  }
-        else
 	  {
 	    for(j=i;j<fIntervalNumber;j++)
@@ -529,5 +527,6 @@
    //  result *= (1-exp(-be2/betaBohr2)) ;
    result *= (1-exp(-be4/betaBohr4)) ;
-   if(fDensity >= 0.1)
+   //   if(fDensity >= 0.1)
+   if(x8 > 0.)
    { 
       result /= x8 ;
@@ -586,10 +585,12 @@
    dNdxC *= (1-exp(-be4/betaBohr4)) ;
 
-   if(fDensity >= 0.1)
-   { 
-      modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 
+   //   if(fDensity >= 0.1)
+   // { 
+   modul2 = (1.0 + fRePartDielectricConst[i])*(1.0 + fRePartDielectricConst[i]) + 
                     fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
-      dNdxC /= modul2 ;
-   }
+   if(modul2 > 0.)
+     {
+       dNdxC /= modul2 ;
+     }
    return dNdxC ;
 
@@ -626,10 +627,12 @@
    dNdxP *= (1-exp(-be4/betaBohr4)) ;
 
-   if( fDensity >= 0.1 )
-   { 
-     modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
-        fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
-     dNdxP /= modul2 ;
-   }
+//   if( fDensity >= 0.1 )
+//   { 
+   modul2 = (1 + fRePartDielectricConst[i])*(1 + fRePartDielectricConst[i]) + 
+     fImPartDielectricConst[i]*fImPartDielectricConst[i] ;
+   if(modul2 > 0.)
+     { 
+       dNdxP /= modul2 ;
+     }
    return dNdxP ;
 
@@ -800,5 +803,5 @@
 G4double G4PAIySection::SumOverInterCerenkov( G4int i )
 {         
-   G4double x0,x1,y0,yy1,a,b,c,result ;
+   G4double x0,x1,y0,yy1,a,c,result ;
 
    x0  = fSplineEnergy[i] ;
@@ -811,5 +814,6 @@
    c = x1/x0;
    a = log10(yy1/y0)/log10(c) ;
-   b = y0/pow(x0,a) ;
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
 
    a += 1.0 ;
@@ -833,5 +837,5 @@
 G4double G4PAIySection::SumOverInterPlasmon( G4int i )
 {         
-   G4double x0,x1,y0,yy1,a,b,c,result ;
+   G4double x0,x1,y0,yy1,a,c,result ;
 
    x0  = fSplineEnergy[i] ;
@@ -841,6 +845,7 @@
    c =x1/x0;
    a = log10(yy1/y0)/log10(c) ;
-   // b = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a) ;
+
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
 
    a += 1.0 ;
@@ -864,5 +869,5 @@
                                        G4double en0    )
 {               
-   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
+   G4double x0,x1,y0,yy1,a,c,d,e0,result ;
 
    e0 = en0 ;
@@ -875,6 +880,7 @@
    d = e0/x0;   
    a = log10(yy1/y0)/log10(x1/x0) ;
-   // b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a);  // pow(10.,b) ;
+
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
    
    a += 1 ;
@@ -933,5 +939,5 @@
                                        G4double en0    )
 {               
-   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
+   G4double x0,x1,y0,yy1,a,c,d,e0,result ;
 
    e0 = en0 ;
@@ -944,6 +950,7 @@
    d = e0/x0;   
    a = log10(yy1/y0)/log10(x1/x0) ;
-   // b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a);  // pow(10.,b) ;
+   
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
    
    a += 2 ;
@@ -964,6 +971,7 @@
    d = e0/x0;   
    a = log10(yy1/y0)/log10(x1/x0) ;
-   //  b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a) ;
+
+   if(a < 20.) b = y0/pow(x0,a) ;
+
    a += 2 ;
    if(a == 0) 
@@ -987,5 +995,5 @@
                                              G4double en0    )
 {               
-   G4double x0,x1,y0,yy1,a,b,e0,c,d,result ;
+   G4double x0,x1,y0,yy1,a,e0,c,d,result ;
 
    e0 = en0 ;
@@ -996,11 +1004,12 @@
 
    //  G4cout<<G4endl;
-   //  G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
+   //G4cout<<"SumBordC, i = "<<i<<"; en0 = "<<en0<<"; x0 ="<<x0<<"; x1 = "<<x1
    //     <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
    c = x1/x0 ;
    d = e0/x0 ;
    a = log10(yy1/y0)/log10(c) ;
-   // b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a); // pow(10.,b0) ;   
+ 
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
    
    a += 1.0 ;
@@ -1012,5 +1021,5 @@
    else         fIntegralCerenkov[0] += y0*(x0*x0 - e0*e0*pow(d,a-2))/a ;
 
-// G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "<<b<<"; result = "<<result<<G4endl;
+   //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;
    
    x0  = fSplineEnergy[i - 1] ;
@@ -1019,5 +1028,5 @@
    yy1 = fdNdxCerenkov[i - 2] ;
 
-   // G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
+   //G4cout<<"x0 ="<<x0<<"; x1 = "<<x1
    //    <<"; y0 = "<<y0<<"; yy1 = "<<yy1<<G4endl;
 
@@ -1025,6 +1034,12 @@
    d = e0/x0 ;
    a  = log10(yy1/y0)/log10(x1/x0) ;
-   // b0 = log10(y0) - a*log10(x0) ;
-   b  =  y0/pow(x0,a);  // pow(10.,b0) ;
+  
+   //   G4cout << "a= " << a << G4endl;
+   if(a < 20.) b = y0/pow(x0,a) ;
+
+   if(a > 20.0) b = 0.0;
+   else         b = y0/pow(x0,a);  // pow(10.,b0) ;
+
+   //G4cout << "b= " << b << G4endl;
 
    a += 1.0 ;
@@ -1032,10 +1047,10 @@
    else         result += y0*(e0*pow(d,a-1) - x0 )/a ;
    a += 1.0 ;
+   //G4cout << "result= " << result << G4endl;
 
    if( a == 0 )   fIntegralCerenkov[0] += b*log(e0/x0) ;
    else           fIntegralCerenkov[0] += y0*(e0*e0*pow(d,a-2) - x0*x0)/a ;
 
-   // G4cout<<"a = "<<a<<"; b0 = "<<b0<<"; b = "
-   // <<b<<"; result = "<<result<<G4endl;    
+   //G4cout<<"a = "<<a<<"; b = "<<b<<"; result = "<<result<<G4endl;    
 
    return result ;
@@ -1051,5 +1066,5 @@
                                              G4double en0    )
 {               
-   G4double x0,x1,y0,yy1,a,b,c,d,e0,result ;
+   G4double x0,x1,y0,yy1,a,c,d,e0,result ;
 
    e0 = en0 ;
@@ -1062,6 +1077,7 @@
    d = e0/x0 ;   
    a = log10(yy1/y0)/log10(c) ;
-   //  b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a); //pow(10.,b) ;
+
+   G4double b = 0.0;
+   if(a < 20.) b = y0/pow(x0,a) ;
    
    a += 1.0 ;
@@ -1081,6 +1097,6 @@
    d = e0/x0 ;
    a = log10(yy1/y0)/log10(c) ;
-   // b0 = log10(y0) - a*log10(x0) ;
-   b = y0/pow(x0,a);// pow(10.,b0) ;
+ 
+   if(a < 20.) b = y0/pow(x0,a) ;
 
    a += 1.0 ;
Index: trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4PEEffectModel.cc,v 1.8 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4PSTARStopping.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PSTARStopping.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PSTARStopping.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4PSTARStopping.cc,v 1.8 2008/11/24 18:28:09 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 
 //---------------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4PairProductionRelModel.cc,v 1.3 2009/05/15 17:12:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4PhotoElectricEffect.cc,v 1.42 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
Index: trunk/source/processes/electromagnetic/standard/src/G4PolarizedComptonScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PolarizedComptonScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PolarizedComptonScattering.cc	(revision 1196)
@@ -26,5 +26,5 @@
 //
 // $Id: G4PolarizedComptonScattering.cc,v 1.18 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 //
Index: trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4UniversalFluctuation.cc,v 1.22 2009/03/20 18:11:23 urban Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc	(revision 1196)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4UrbanMscModel.cc,v 1.90 2009/04/29 13:30:22 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4UrbanMscModel.cc,v 1.91 2009/07/20 18:41:34 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -122,5 +122,5 @@
   lambdalimit   = 1.*mm;
   fr            = 0.02;
-  facsafety     = 0.3;
+  //  facsafety     = 0.3;
   taubig        = 8.0;
   tausmall      = 1.e-16;
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc	(revision 1196)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4UrbanMscModel2.cc,v 1.26 2009/05/19 06:26:10 urban Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4UrbanMscModel2.cc,v 1.27 2009/07/20 18:41:34 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -119,5 +119,5 @@
   lambdalimit   = 1.*mm;
   fr            = 0.02;
-  facsafety     = 0.3;
+  //facsafety     = 0.3;
   taubig        = 8.0;
   tausmall      = 1.e-16;
@@ -462,7 +462,8 @@
   presafety = sp->GetSafety();
 
-  //    G4cout << "G4UrbanMscModel2::ComputeTruePathLengthLimit tPathLength= " 
-  //	   <<tPathLength<<" safety= " << presafety
-  //	   << " range= " <<currentRange<<G4endl;
+  //G4cout << "G4Urban2::StepLimit tPathLength= " 
+  //	 <<tPathLength<<" safety= " << presafety
+  //       << " range= " <<currentRange<< " lambda= "<<lambda0
+  //	 << " Alg: " << steppingAlgorithm <<G4endl;
 
   // far from geometry boundary
@@ -506,5 +507,6 @@
           tlimitmin = 10.*stepmin;
           if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
-
+	  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
+	  //	 << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
           // constraint from the geometry
           if((geomlimit < geombig) && (geomlimit > geommin))
@@ -535,6 +537,6 @@
       if(tlimit > tgeom) tlimit = tgeom;
 
-      //  G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
-      //     << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
+      //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
+      //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
 
       // shortcut
@@ -627,5 +629,5 @@
 	}
     }
-  // G4cout << "tPathLength= " << tPathLength << "  geomlimit= " << geomlimit 
+  //G4cout << "tPathLength= " << tPathLength 
   //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
   return tPathLength ;
@@ -713,5 +715,5 @@
 
   if(zPathLength > lambda0) zPathLength = lambda0;
-
+  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
   return zPathLength;
 }
@@ -744,4 +746,5 @@
   }
   if(tPathLength < geomStepLength) tPathLength = geomStepLength;
+  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
 
   return tPathLength;
@@ -799,5 +802,5 @@
 
     G4double r = SampleDisplacement();
-    /*
+    /*    
     G4cout << "G4UrbanMscModel2::SampleSecondaries: e(MeV)= " << kineticEnergy
 	   << " sinTheta= " << sth << " r(mm)= " << r
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4UrbanMscModel90.cc,v 1.13 2009/04/10 18:10:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WaterStopping.cc,v 1.16 2009/05/15 19:04:21 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4WaterStopping.cc,v 1.18 2009/06/19 10:39:48 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 
 //---------------------------------------------------------------------------
@@ -67,5 +67,5 @@
 {
   G4double res = 0.0;
-  if(iz < 3 || iz > 18) return res; 
+  if((iz > 26) || (iz < 3) || (iz > 18 && iz < 26)) return res; 
   G4bool b;
   G4int idx = iz - 3;
@@ -99,8 +99,8 @@
   G4int i;
   //..List of ions
-  G4int zz[16] = {3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18};
-  G4int aa[16] = {7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28,31,32, 35,40};
-  // G4double A_Ion[17] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948};
-  for(i=0; i<16; i++) {
+  G4int zz[17] = {3, 4, 5, 6, 7, 8, 9, 10,11,12,13,14,15,16,17,18,26};
+  G4int aa[17] = {7, 9, 11, 12, 14, 16, 19, 20, 23, 24, 27, 28,31,32, 35,40,56};
+  // G4double A_Ion[17] = {6.941,9.0122,10.811,12.011,14.007,15.999,18.998,20.180,22.990,24.305,26.982,28.086,30.974,32.065,35.453,39.948,55.845};
+  for(i=0; i<17; i++) {
     Z[i] = zz[i];
     A[i] = G4double(aa[i]);
@@ -145,6 +145,9 @@
   AddData(E,G4_WATER_Ar,factor);
 
+  G4double G4_WATER_Fe [53]={6.5394, 7.3060, 8.7367, 10.0690, 11.3310, 12.5470, 13.7280, 14.8780, 15.9980, 21.1160, 25.4850, 29.1260, 32.0640, 36.1770, 38.6920, 40.2760, 41.2950, 41.9530, 42.3710, 42.6210, 42.5910, 41.6900, 40.5190, 39.2690, 36.8000, 34.5040, 32.4190, 30.5410, 28.8480, 27.3170, 25.9310, 20.6170, 17.0680, 14.5540, 12.6930, 10.1410, 8.4892, 7.3402, 6.4976, 5.8545, 5.3479, 4.9387, 3.6892, 3.0503, 2.6620, 2.4014, 2.0756, 1.8825, 1.7569, 1.6702, 1.6079, 1.5619, 1.5267};
+  AddData(E,G4_WATER_Fe,factor);
+
   if(corr) {
-    for(i=0; i<16; i++) {corr->AddStoppingData(Z[i], aa[i], "G4_WATER", dedx[i]);}
+    for(i=0; i<17; i++) {corr->AddStoppingData(Z[i], aa[i], "G4_WATER", dedx[i]);}
   }
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WentzelVIModel.cc,v 1.32 2009/05/10 16:09:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4WentzelVIModel.cc,v 1.37 2009/10/28 10:14:13 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -135,5 +135,5 @@
   // reset parameters
   SetupParticle(p);
-  tkin = targetZ = mom2 = DBL_MIN;
+  tkin = targetZ = mom2 = 0.0;
   ecut = etag = DBL_MAX;
   currentRange = 0.0;
@@ -158,9 +158,9 @@
 {
   SetupParticle(p);
-  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
-  SetupKinematic(ekin, cutEnergy);
-  SetupTarget(Z, ekin);
+  if(kinEnergy < lowEnergyLimit) return 0.0;
+  SetupKinematic(kinEnergy, cutEnergy);
+  SetupTarget(Z, kinEnergy);
   G4double xsec = ComputeTransportXSectionPerAtom();
-  /*  
+  /*   
   G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2 
 	 << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn
@@ -392,5 +392,7 @@
   //	 << particle->GetParticleName() << G4endl;
   G4double kinEnergy = dynParticle->GetKineticEnergy();
-  if(kinEnergy <= DBL_MIN || tPathLength <= DBL_MIN) return;
+
+  // ignore scattering for zero step length and enegy below the limit
+  if(kinEnergy < lowEnergyLimit || tPathLength <= DBL_MIN) return;
 
   G4double ekin = preKinEnergy;
@@ -423,5 +425,4 @@
 
     // for low-energy e-,e+ no limit
-    ekin = std::max(ekin, lowEnergyLimit);
     SetupKinematic(ekin, cut);
   
@@ -431,6 +432,6 @@
       xsec = ComputeXSectionPerVolume();
 
-      if(xtsec > DBL_MIN) x1 = 0.5*tPathLength*xtsec;
-      else                x1 = 0.0;
+      if(xtsec > 0.0) x1 = 0.5*tPathLength*xtsec;
+      else            x1 = 0.0;
 
       /*      
@@ -452,6 +453,6 @@
   // cost is sampled ------------------------------
   G4double cost = 1.0 - 2.0*z;
-  if(cost < -1.0) cost = -1.0;
-  else if(cost > 1.0) cost = 1.0;
+  //  if(cost < -1.0) cost = -1.0;
+  // else if(cost > 1.0) cost = 1.0;
   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
 
@@ -472,5 +473,5 @@
 
   // sample MSC scattering for large angle
-  // extra central scattering for holf step
+  // extra central scattering for half step
   if(largeAng) {
     isscat = true;
@@ -480,5 +481,5 @@
     } while (z > 1.0); 
     cost = 1.0 - 2.0*z;
-    if(std::abs(cost) > 1.0) cost = 1.0;
+    //if(std::fabs(cost) > 1.0) cost = 1.0;
 
     sint = sqrt((1.0 - cost)*(1.0 + cost));
@@ -491,5 +492,5 @@
   }
 
-  // sample Reserford scattering for large angle
+  // sample Rutherford scattering for large angle
   if(xsec > DBL_MIN) {
     G4double t = tPathLength;
@@ -536,5 +537,5 @@
         if(zz1 < 1.0) {
 	  isscat = true;
-	  //G4cout << "Reserford zz1= " << zz1 << " t= " << t << G4endl;
+	  //G4cout << "Rutherford zz1= " << zz1 << " t= " << t << G4endl;
 	  sint = sqrt((1.0 - zz1)*(1.0 + zz1));
 	  //G4cout << "sint= " << sint << G4endl;
@@ -637,5 +638,5 @@
     if(cosThetaMin > cosnm) {
 
-      // Reserford part
+      // Rutherford part
       G4double s  = screenZ*formfactA;
       G4double z1 = 1.0 - cosnm + screenZ;
@@ -711,5 +712,5 @@
       zcorr  += 0.5*x*s;
 
-      // Reserford part
+      // Rutherford part
       G4double z1 = 1.0 - cosTetMaxElec + screenZ;
       G4double z2 = (cosThetaMin - cosTetMaxElec)/x1; 
@@ -730,5 +731,5 @@
       zcorr  += 0.5*x*s;
 
-      // Reserford part
+      // Rutherford part
       s  = screenZ*formfactA;
       G4double w  = 1.0 + 2.0*s;
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eBremsstrahlung.cc,v 1.56 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eBremsstrahlungModel.cc,v 1.44 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eBremsstrahlungRelModel.cc,v 1.14 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eCoulombScatteringModel.cc,v 1.69 2009/05/10 16:09:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// $Id: G4eCoulombScatteringModel.cc,v 1.78 2009/10/28 10:14:13 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -39,4 +39,5 @@
 //
 // Modifications:
+//
 // 01.08.06 V.Ivanchenko extend upper limit of table to TeV and review the
 //          logic of building - only elements from G4ElementTable
@@ -45,4 +46,6 @@
 // 09.10.07 V.Ivanchenko reorganized methods, add cut dependence in scattering off e- 
 // 09.06.08 V.Ivanchenko add SelectIsotope and sampling of the recoil ion 
+// 16.06.09 C.Consolandi fixed computation of effective mass
+//
 //
 // Class Description:
@@ -63,4 +66,6 @@
 #include "G4Proton.hh"
 #include "G4ParticleTable.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4NucleiProperties.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -87,10 +92,10 @@
   currentMaterial = 0; 
   currentElement  = 0;
-  lowEnergyLimit = keV;
+  lowEnergyLimit  = 0.1*keV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff  = twopi*p0*p0;
   tkin = targetZ = mom2 = DBL_MIN;
   elecXSection = nucXSection = 0.0;
-  recoilThreshold = 100.*keV;
+  recoilThreshold = 0.*keV;
   ecut = DBL_MAX;
   particle = 0;
@@ -130,8 +135,9 @@
   ecut = etag = DBL_MAX;
   cosThetaMin = cos(PolarAngleLimit());
-  currentCuts = &cuts;
+  pCuts = G4ProductionCutsTable::GetProductionCutsTable()->GetEnergyCutsVector(3);
   //G4cout << "!!! G4eCoulombScatteringModel::Initialise for " 
   //	 << p->GetParticleName() << "  cos(TetMin)= " << cosThetaMin 
   //	 << "  cos(TetMax)= " << cosThetaMax <<G4endl;
+  // G4cout << "cut0= " << cuts[0] << "  cut1= " << cuts[1] << G4endl;
   if(!isInitialised) {
     isInitialised = true;
@@ -185,8 +191,8 @@
   G4double xsec = 0.0;
   SetupParticle(p);
-  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
-  SetupKinematic(ekin, cutEnergy);
+  if(kinEnergy < lowEnergyLimit) return xsec;
+  SetupKinematic(kinEnergy, cutEnergy);
   if(cosTetMaxNuc < cosTetMinNuc) {
-    SetupTarget(Z, ekin);
+    SetupTarget(Z, kinEnergy);
     xsec = CrossSectionPerAtom();  
   }
@@ -196,6 +202,5 @@
 	 << " cosTetMaxElec= " << cosTetMaxElec
 	 << " screenZ= " << screenZ
-	 << " formfactA= " << formfactA
-	 << " cosTetMaxHad= " << cosTetMaxHad << G4endl;
+	 << " formfactA= " << formfactA << G4endl;
   */
   return xsec;  
@@ -207,5 +212,9 @@
 {
   // This method needs initialisation before be called
-  G4double fac = coeff*targetZ*chargeSquare*kinFactor;
+  //G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2;
+
+  G4double meff = targetMass/(mass+targetMass);
+  G4double fac  = coeff*targetZ*chargeSquare*invbeta2/(mom2*meff*meff);
+
   elecXSection = 0.0;
   nucXSection  = 0.0;
@@ -257,16 +266,21 @@
 {
   G4double kinEnergy = dp->GetKineticEnergy();
-  if(kinEnergy <= DBL_MIN) return;
+  if(kinEnergy < lowEnergyLimit) return;
   DefineMaterial(couple);
   SetupParticle(dp->GetDefinition());
-  G4double ekin = std::max(lowEnergyLimit, kinEnergy);
-  SetupKinematic(ekin, cutEnergy);
+
+  SetupKinematic(kinEnergy, cutEnergy);
   //G4cout << "G4eCoulombScatteringModel::SampleSecondaries e(MeV)= " 
-  //	 << kinEnergy << "  " << particle->GetParticleName() << G4endl;
+  //	 << kinEnergy << "  " << particle->GetParticleName() 
+  //	 << " cut= " << cutEnergy<< G4endl;
  
   // Choose nucleus
-  currentElement = SelectRandomAtom(couple,particle,ekin,cutEnergy,ekin);
-
-  SetupTarget(currentElement->GetZ(),ekin);
+  currentElement = SelectRandomAtom(couple,particle,
+				    kinEnergy,cutEnergy,kinEnergy);
+
+  SetupTarget(currentElement->GetZ(),kinEnergy);
+
+  G4int ia = SelectIsotopeNumber(currentElement);
+  targetMass = G4NucleiProperties::GetNuclearMass(ia, iz);
   
   G4double cost = SampleCosineTheta();
@@ -276,5 +290,5 @@
   G4double sint = sqrt(z1*(1.0 + cost));
   
-  //G4cout<<"## Sampled sint= " << sint << "  Z= " << targetZ 
+  //G4cout<<"## Sampled sint= " << sint << "  Z= " << targetZ << " A= " << ia
   //	<< "  screenZ= " << screenZ << " cn= " << formfactA << G4endl;
   
@@ -289,18 +303,26 @@
   // recoil sampling assuming a small recoil
   // and first order correction to primary 4-momentum
-  if(lowEnergyLimit < kinEnergy) {
-    G4int ia = SelectIsotopeNumber(currentElement);
-    G4double Trec = z1*mom2/(amu_c2*G4double(ia));
-
-    if(Trec > recoilThreshold) {
-      G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
-      Trec = z1*mom2/ion->GetPDGMass();
-      if(Trec < kinEnergy) {
-	G4ThreeVector dir = (direction - newDirection).unit();
-	G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, Trec);
-	fvect->push_back(newdp);
-	fParticleChange->SetProposedKineticEnergy(kinEnergy - Trec);
-      }
-    }
+  G4double q2   = 2*z1*mom2;
+  G4double trec = q2/(sqrt(targetMass*targetMass + q2) + targetMass);
+  G4double finalT = kinEnergy - trec; 
+  //G4cout<<"G4eCoulombScatteringModel: finalT= "<<finalT<<" Trec= "<<trec<<G4endl;
+  if(finalT <= lowEnergyLimit) { 
+    trec = kinEnergy;  
+    finalT = 0.0;
+  } 
+
+  fParticleChange->SetProposedKineticEnergy(finalT);
+  G4double tcut = recoilThreshold;
+  if(pCuts) { tcut= std::max(tcut,(*pCuts)[currentMaterialIndex]); }
+
+  if(trec > tcut) {
+    G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
+    G4ThreeVector dir = (direction*sqrt(mom2) - 
+			 newDirection*sqrt(finalT*(2*mass + finalT))).unit();
+    G4DynamicParticle* newdp = new G4DynamicParticle(ion, dir, trec);
+    fvect->push_back(newdp);
+  } else {
+    fParticleChange->ProposeLocalEnergyDeposit(trec);
+    fParticleChange->ProposeNonIonizingEnergyDeposit(trec);
   }
  
@@ -324,12 +346,12 @@
   }
 
-  /*
+  /*  
   G4cout << "SampleCost: e(MeV)= " << tkin 
-  	 << " ctmin= " << cosThetaMin
-  	 << " ctmaxN= " << cosTetMaxNuc
-  	 << " ctmax= " << costm
-  	 << " Z= " << targetZ << " A= " << targetA
+  	 << " 1-ctmaxN= " << 1. - cosTetMinNuc
+  	 << " 1-ctmax= " << 1. - costm
+  	 << " Z= " << targetZ 
   	 << G4endl;
   */
+
   if(costm >= cosTetMinNuc) return 2.0; 
 
@@ -343,6 +365,11 @@
   } while ( G4UniformRand() > grej*grej );  
 
-  //G4cout << "z= " << z1 << " cross= " << nucXSection/barn 
-  // << " crossE= " << elecXSection/barn << G4endl;
+  if(mass > MeV) {
+    if(G4UniformRand() > (1. - z1*0.5)/(1.0 + z1*sqrt(mom2)/targetMass)) {
+      return 2.0;
+    }
+  }
+  //G4cout << "z1= " << z1 << " cross= " << nucXSection/barn 
+  //	 << " crossE= " << elecXSection/barn << G4endl;
 
   return 1.0 - z1;
Index: trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eIonisation.cc,v 1.57 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4eMultipleScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eMultipleScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eMultipleScattering.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eMultipleScattering.cc,v 1.7 2008/10/23 17:55:20 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eMultipleScattering.cc,v 1.10 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -----------------------------------------------------------------------------
@@ -45,5 +45,5 @@
 
 #include "G4eMultipleScattering.hh"
-#include "G4UrbanMscModel2.hh"
+#include "G4UrbanMscModel92.hh"
 #include "G4MscStepLimitType.hh"
 #include "G4Electron.hh"
@@ -58,5 +58,4 @@
 {
   isInitialized = false;  
-  SetRangeFactor(0.04);
 }
 
@@ -76,34 +75,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4eMultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)
+void G4eMultipleScattering::InitialiseProcess(const G4ParticleDefinition*)
 {
-  // Modification of parameters between runs
-  if(isInitialized) {
-    if (p->GetParticleType() != "nucleus") {
-      mscUrban->SetStepLimitType(StepLimitType());
-      mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-      mscUrban->SetSkin(Skin());
-      mscUrban->SetRangeFactor(RangeFactor());
-      mscUrban->SetGeomFactor(GeomFactor());
-    }
-    return;
-  }
-
-  // defaults for ions, which cannot be overwritten
-  if (p->GetParticleType() == "nucleus") {
-    SetStepLimitType(fMinimal);
-    SetLateralDisplasmentFlag(false);
-    SetBuildLambdaTable(false);
-  }
+  if(isInitialized) return;
 
   // initialisation of parameters - defaults for particles other
   // than ions can be overwritten by users
-  mscUrban = new G4UrbanMscModel2();
-  mscUrban->SetStepLimitType(StepLimitType());
-  mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-  mscUrban->SetSkin(Skin());
-  mscUrban->SetRangeFactor(RangeFactor());
-  mscUrban->SetGeomFactor(GeomFactor());
-
+  G4VMscModel* mscUrban = new G4UrbanMscModel92();
   AddEmModel(1,mscUrban);
   isInitialized = true;
Index: trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eeToTwoGammaModel.cc,v 1.15 2009/04/09 18:41:18 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4eplusAnnihilation.cc,v 1.30 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4hIonisation.cc,v 1.82 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
Index: trunk/source/processes/electromagnetic/standard/src/G4hMultipleScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4hMultipleScattering.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4hMultipleScattering.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4hMultipleScattering.cc,v 1.13 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4hMultipleScattering.cc,v 1.16 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -----------------------------------------------------------------------------
@@ -59,5 +59,4 @@
 {
   isInitialized = false;  
-  isIon         = false;
   SetStepLimitType(fMinimal);
 }
@@ -77,36 +76,8 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4hMultipleScattering::InitialiseProcess(const G4ParticleDefinition* p)
+void G4hMultipleScattering::InitialiseProcess(const G4ParticleDefinition*)
 {
-  // Modification of parameters between runs
-  if(isInitialized) {
-    if (p->GetParticleType() != "nucleus" && p->GetPDGMass() < GeV) {
-      mscUrban->SetStepLimitType(StepLimitType());
-      mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-      mscUrban->SetSkin(Skin());
-      mscUrban->SetRangeFactor(RangeFactor());
-      mscUrban->SetGeomFactor(GeomFactor());
-    }
-    return;
-  }
-
-  // defaults for ions, which cannot be overwritten
-  if (p->GetParticleType() == "nucleus" || p->GetPDGMass() > GeV) {
-    SetStepLimitType(fMinimal);
-    SetLateralDisplasmentFlag(false);
-    SetBuildLambdaTable(false);
-    if(p->GetParticleType() == "nucleus") isIon = true;
-  }
-
-  // initialisation of parameters
-  G4String part_name = p->GetParticleName();
-  mscUrban = new G4UrbanMscModel90();
-
-  mscUrban->SetStepLimitType(StepLimitType());
-  mscUrban->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-  mscUrban->SetSkin(Skin());
-  mscUrban->SetRangeFactor(RangeFactor());
-  mscUrban->SetGeomFactor(GeomFactor());
-
+  if(isInitialized) return;
+  G4VMscModel* mscUrban = new G4UrbanMscModel90();
   AddEmModel(1,mscUrban);
   isInitialized = true;
@@ -118,8 +89,8 @@
 {
   G4cout << "      RangeFactor= " << RangeFactor()
-	 << ", step limit type: " << StepLimitType()
-         << ", lateralDisplacement: " << LateralDisplasmentFlag()
-	 << ", skin= " << Skin()  
-    //	 << ", geomFactor= " << GeomFactor()  
+	 << ", stepLimitType: " << StepLimitType()
+         << ", latDisplacement: " << LateralDisplasmentFlag()
+	 << ", skin= " << Skin()
+   	 << ", geomFactor= " << GeomFactor()  
 	 << G4endl;
 }
@@ -127,27 +98,2 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4hMultipleScattering::AlongStepGetPhysicalInteractionLength(
-                             const G4Track& track,
-                             G4double,
-                             G4double currentMinimalStep,
-                             G4double& currentSafety,
-                             G4GPILSelection* selection)
-{
-  // get Step limit proposed by the process
-  valueGPILSelectionMSC = NotCandidateForSelection;
-
-  G4double escaled = track.GetKineticEnergy();
-  if(isIon) escaled *= track.GetDynamicParticle()->GetMass()/proton_mass_c2;
-
-  G4double steplength = GetMscContinuousStepLimit(track,
-						  escaled,
-						  currentMinimalStep,
-						  currentSafety);
-  // G4cout << "StepLimit= " << steplength << G4endl;
-  // set return value for G4GPILSelection
-  *selection = valueGPILSelectionMSC;
-  return  steplength;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
Index: trunk/source/processes/electromagnetic/standard/src/G4ionGasIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ionGasIonisation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4ionGasIonisation.cc	(revision 1196)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4ionGasIonisation.cc,v 1.14 2008/09/12 16:26:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ionGasIonisation.cc,v 1.15 2009/11/10 11:50:30 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -58,4 +58,8 @@
   : G4ionIonisation(name)
 {
+  G4cout << G4endl;
+  G4cout << "!!! G4ionGasIonisation class is obsolete and may be removed for the next major Geant4 release !!!" << G4endl;
+  G4cout << "!!! Please use G4ionIonisation for ions !!!" << G4endl;
+  G4cout << G4endl;
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc	(revision 1196)
@@ -25,5 +25,5 @@
 //
 // $Id: G4ionIonisation.cc,v 1.66 2009/02/20 12:06:37 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
 //
 // -------------------------------------------------------------------
