Index: trunk/source/processes/electromagnetic/standard/History
===================================================================
--- trunk/source/processes/electromagnetic/standard/History	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/History	(revision 1055)
@@ -1,3 +1,3 @@
-$Id: History,v 1.430 2008/11/24 18:28:09 vnivanch Exp $
+$Id: History,v 1.452 2009/05/18 12:31:46 vnivanch Exp $
 -------------------------------------------------------------------
 
@@ -17,4 +17,80 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+
+18 May 09:  V.Ivant (emstand-V09-02-12)
+- G4UrbanMscModel2  - (L.Urban) correction in tail tuning using MuScat data
+
+15 May 09:  V.Ivant (emstand-V09-02-11)
+- G4WaterStopping  - fixed Ar data
+- G4PairProductionRelModel - (A.Schaelicke) new relativistic model for 
+                              gamma conversion
+- G4UrbanMscModel2  - (L.Urban) new version of theta0 and tail tuning
+
+10 May 09: V.Ivant (emstand-V09-02-10)
+- G4WentzelVIModel, G4eCoulombScattering, G4CoulombScattering - added relativistic
+                      factor to Reserford cross section; set default limit on
+                      kinetic energy of the recoil 100 keV
+
+28 April 09: V.Ivant (emstand-V09-02-09)
+- G4UrbanMscModel2  - (L.Urban) new tuning for the central part and for the 
+                      tail of the angular distribution using the old e- 
+                      scattering data only (Phys. Rev. 84 (1951) 634; 
+                      Phys. Rev. 61 (1942) 254) 
+                    - corrected logic in ComputeTruePathLengthLimit method 
+                      for type=fUseDistanceToBoundary
+- G4UrbanMscModel   - frozen version of G4UrbanMscModel2 of g4 9.2
+- G4WentzelVIModel, G4eCoulombScattering  - reduce low-limit from 1 keV to 0.1 keV 
+                      to provide smooth transport cross section table
+
+23 April 09: V.Ivant (emstand-V09-02-08)
+- G4BetheBlochModel - do not use pointer to GenericIon introduced in the
+                      previous tag due to problem of simple PhysLists without ions
+
+20 April 09: V.Ivant (emstand-V09-02-07)
+- G4BetheBlochModel - fixed and simplified initialisation for ions needed for
+                      the new G4IonParametrisedLossModel of low-energy package
+- G4GoudsmitSaundersonMscModel - (O.Kadri) cleanup: discarded no scattering and 
+                      single scattering theta sampling from SampleCosineTheta() 
+                      which means the splitting step into two sub-steps occur 
+                      only for msc regime  
+
+12 April 09: V.Ivant (emstand-V09-02-06)
+- Simplified initialisation of all models
+- G4UrbanMscModel, G4UrbanMscModel2, G4UrbanMscModel90 - use methods
+   of G4VMscModel for interface to geometry
+
+07 April 09: V.Ivant (emstand-V09-02-05)
+- G4IonFluctuations - T.Toshito removed extra phenomenological factor
+                      in fluctuation width
+- G4HeatedKleinNishinaCompton - V.Grichine added a prototype model for
+                                plasma 
+
+21 March 09: V.Ivant (emstand-V09-02-04)
+- G4UniversalFluctuation - L.Urban introduce modification in width 
+    correction, the dependence of the correction on energy deposition 
+    at previous steps is removed (addresses T2K report)
+
+20 March 09: V.Ivant (emstand-V09-02-03)
+- G4GoudsmitSaundersonMscModel fixed compillation problem
+
+12 March 09: V.Ivant (emstand-V09-02-02)
+- G4GoudsmitSaundersonMscModel fixed compillation problem
+- G4UniversalFluctuation - add temporary fix for T2K report
+
+05 March 09: V.Ivant (emstand-V09-02-01)
+- New G4GoudsmitSaundersonMscModel is added
+- G4WentzelVIModel, G4eCoulombScattringModel:
+    o substitute scaling of low-energy limit by setting 1 keV for 
+      all particles;
+    o use EGSnrc form of screening parameter (second order correction) 
+
+20 February 09: V.Ivant (emstand-V09-02-00)
+- Move all virtual methods from inline to source
+G4PEEffectModel - substitute ComputeMeanFreePath by CrossSectionPerVolume 
+                  (minor CPU speadup for compound materials)
+G4PAIModel, G4PAIPhotonModel - remove usage of random numbers at 
+                   initialisation (potential non-reproducibility)
+G4WentzelVIModel - use generic methods of G4VMscModel to access safety
+                   and other geometry information 
 
 24 November 08: V.Ivant (emstand-V09-01-45)
Index: trunk/source/processes/electromagnetic/standard/include/G4BetheBlochModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4BetheBlochModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4BetheBlochModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BetheBlochModel.hh,v 1.16 2008/10/22 16:00:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BetheBlochModel.hh,v 1.20 2009/04/23 17:44:43 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -64,9 +64,8 @@
 
 #include "G4VEmModel.hh"
+#include "G4NistManager.hh"
 
 class G4EmCorrections;
 class G4ParticleChangeForLoss;
-class G4NistManager;
-
 
 class G4BetheBlochModel : public G4VEmModel
@@ -82,6 +81,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*,
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
+				const G4MaterialCutsCouple*);
 			
   virtual G4double ComputeCrossSectionPerElectron(
@@ -131,10 +130,14 @@
 protected:
 
-  G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-			      G4double kinEnergy);
+  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
+				      G4double kinEnergy);
 
 private:
 
-  void SetParticle(const G4ParticleDefinition* p);
+  inline void SetupParameters();
+
+  inline void SetParticle(const G4ParticleDefinition* p);
+
+  inline void SetGenericIon(const G4ParticleDefinition* p);
 
   // hide assignment operator
@@ -166,16 +169,45 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline G4double G4BetheBlochModel::MaxSecondaryEnergy(
-          const G4ParticleDefinition* pd,
-                G4double kinEnergy) 
-{
-  if(isIon) SetParticle(pd);
-  G4double tau  = kinEnergy/mass;
-  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
-                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
-  return std::min(tmax,tlimit);
+inline void G4BetheBlochModel::SetupParameters()
+{
+  mass = particle->GetPDGMass();
+  spin = particle->GetPDGSpin();
+  G4double q = particle->GetPDGCharge()/eplus;
+  chargeSquare = q*q;
+  ratio = electron_mass_c2/mass;
+  G4double magmom = particle->GetPDGMagneticMoment()*mass/(0.5*eplus*hbar_Planck*c_squared);
+  magMoment2 = magmom*magmom - 1.0;
+  formfact = 0.0;
+  if(particle->GetLeptonNumber() == 0) {
+    G4double x = 0.8426*GeV;
+    if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
+    else if(mass > GeV) {
+      x /= nist->GetZ13(mass/proton_mass_c2);
+      //	tlimit = 51.2*GeV*A13[iz]*A13[iz];
+    }
+    formfact = 2.0*electron_mass_c2/(x*x);
+    tlimit   = 2.0/formfact;
+  }
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+inline void G4BetheBlochModel::SetParticle(const G4ParticleDefinition* p)
+{
+  if(particle != p) {
+    particle = p;
+    if (p->GetPDGCharge()/eplus > 1.5 && p->GetBaryonNumber() > 2) isIon = true;
+    SetupParameters();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline void G4BetheBlochModel::SetGenericIon(const G4ParticleDefinition* p)
+{
+  if(p && particle != p) {
+    if(p->GetParticleName() == "GenericIon") isIon = true;
+  }
+}
+
 #endif
Index: trunk/source/processes/electromagnetic/standard/include/G4BohrFluctuations.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4BohrFluctuations.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4BohrFluctuations.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BohrFluctuations.hh,v 1.3 2007/09/27 13:53:11 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BohrFluctuations.hh,v 1.4 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -79,6 +79,4 @@
   void InitialiseMe(const G4ParticleDefinition*);
 
-protected:
-
 private:
 
@@ -103,25 +101,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-
-inline G4double G4BohrFluctuations::Dispersion(
-                          const G4Material* material,
-                          const G4DynamicParticle* dp,
- 				G4double& tmax,
-			        G4double& length)
-{
-  if(!particle) InitialiseMe(dp->GetDefinition());
-
-  G4double electronDensity = material->GetElectronDensity();
-  kineticEnergy = dp->GetKineticEnergy();
-  G4double etot = kineticEnergy + particleMass;
-  beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
-  G4double siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
-                 * electronDensity * chargeSquare;
-
-  return siga;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 #endif
 
Index: trunk/source/processes/electromagnetic/standard/include/G4BraggIonModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4BraggIonModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4BraggIonModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggIonModel.hh,v 1.11 2008/10/22 16:00:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BraggIonModel.hh,v 1.12 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -77,6 +77,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*,
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
+				const G4MaterialCutsCouple*);
 			
   virtual G4double ComputeCrossSectionPerElectron(
@@ -128,6 +128,6 @@
 protected:
 
-  G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-			      G4double kinEnergy);
+  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
+				      G4double kinEnergy);
 
 private:
@@ -180,17 +180,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-inline G4double G4BraggIonModel::MaxSecondaryEnergy(
-          const G4ParticleDefinition* pd,
-                G4double kinEnergy)
-{
-  if(pd != particle) SetParticle(pd);
-  G4double tau  = kinEnergy/mass;
-  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
-                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
-  return tmax;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
 inline void G4BraggIonModel::SetParticle(const G4ParticleDefinition* p)
 {
Index: trunk/source/processes/electromagnetic/standard/include/G4BraggModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4BraggModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4BraggModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggModel.hh,v 1.12 2008/09/14 17:11:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BraggModel.hh,v 1.13 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -80,6 +80,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*,
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
+				const G4MaterialCutsCouple*);
 			
   virtual G4double ComputeCrossSectionPerElectron(
@@ -131,10 +131,10 @@
 protected:
 
-  G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-			      G4double kinEnergy);
+  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
+				      G4double kinEnergy);
 
 private:
 
-  void SetParticle(const G4ParticleDefinition* p);
+  inline void SetParticle(const G4ParticleDefinition* p);
 
   G4bool HasMaterial(const G4Material* material);
@@ -182,17 +182,4 @@
   G4bool   isInitialised;
 };
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline G4double G4BraggModel::MaxSecondaryEnergy(
-                                            const G4ParticleDefinition* pd,
-                                                  G4double kinEnergy)
-{
-  if(pd != particle) SetParticle(pd);
-  G4double tau  = kinEnergy/mass;
-  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
-                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
-  return tmax;
-}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/include/G4ComptonScattering.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4ComptonScattering.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4ComptonScattering.hh	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4ComptonScattering.hh,v 1.20 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ComptonScattering.hh,v 1.21 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //------------------ G4ComptonScattering physics process -----------------------
@@ -82,9 +82,8 @@
 
   // true for Gamma only.  
-  G4bool IsApplicable(const G4ParticleDefinition&);
+  virtual G4bool IsApplicable(const G4ParticleDefinition&);
   
   // Print few lines of informations about the process: validity range,
   virtual void PrintInfo();
-
 
 protected:
@@ -97,13 +96,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline 
-G4bool G4ComptonScattering::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Gamma::Gamma());
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
   
Index: trunk/source/processes/electromagnetic/standard/include/G4CoulombScattering.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4CoulombScattering.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4CoulombScattering.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScattering.hh,v 1.11 2008/06/13 08:19:43 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4CoulombScattering.hh,v 1.13 2009/05/07 18:41:45 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -59,5 +59,5 @@
 public:
 
-  G4CoulombScattering(const G4String& name = "eCoulombScat");
+  G4CoulombScattering(const G4String& name = "CoulombScat");
 
   virtual ~G4CoulombScattering();
@@ -104,11 +104,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline G4bool G4CoulombScattering::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived());
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 inline void G4CoulombScattering::SetThetaMin(G4double val)
 {
Index: trunk/source/processes/electromagnetic/standard/include/G4GammaConversion.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4GammaConversion.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4GammaConversion.hh	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4GammaConversion.hh,v 1.22 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4GammaConversion.hh,v 1.23 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //
@@ -85,5 +85,5 @@
 
   // true for Gamma only.
-  G4bool IsApplicable(const G4ParticleDefinition&);
+  virtual G4bool IsApplicable(const G4ParticleDefinition&);
 
   // Print few lines of informations about the process: validity range,
@@ -99,12 +99,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4GammaConversion::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Gamma::Gamma());
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
   
Index: trunk/source/processes/electromagnetic/standard/include/G4IonFluctuations.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4IonFluctuations.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4IonFluctuations.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4IonFluctuations.hh,v 1.8 2008/10/22 16:04:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4IonFluctuations.hh,v 1.9 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -84,5 +84,5 @@
 
   // Initialisation prestep
-  inline void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2);
+  void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2);
 
 private:
@@ -117,20 +117,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline
-void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part,
-					     G4double q2)
-{
-  if(part != particle) {
-    particle       = part;
-    particleMass   = part->GetPDGMass();
-    charge         = part->GetPDGCharge()/eplus;
-    chargeSquare   = charge*charge;
-  }
-  effChargeSquare  = q2;
-  uniFluct.SetParticleAndCharge(part, q2);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 #endif
 
Index: trunk/source/processes/electromagnetic/standard/include/G4MollerBhabhaModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4MollerBhabhaModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4MollerBhabhaModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MollerBhabhaModel.hh,v 1.19 2007/05/22 17:34:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4MollerBhabhaModel.hh,v 1.20 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -74,6 +74,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*,
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*,
+				const G4MaterialCutsCouple*);
 				
   virtual G4double ComputeCrossSectionPerElectron(
@@ -109,8 +109,8 @@
 protected:
 
-  G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-			      G4double kinEnergy);
+  virtual G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
+				      G4double kinEnergy);
 			      
-  void SetParticle(const G4ParticleDefinition* p);			      
+  inline void SetParticle(const G4ParticleDefinition* p);			      
 
   const G4ParticleDefinition* particle;
@@ -127,15 +127,15 @@
   G4MollerBhabhaModel & operator=(const  G4MollerBhabhaModel &right);
   G4MollerBhabhaModel(const  G4MollerBhabhaModel&);
+
+  G4bool   isInitialised;
+
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-inline G4double G4MollerBhabhaModel::MaxSecondaryEnergy(
-                                                   const G4ParticleDefinition*,
-                                                         G4double kinEnergy) 
+inline void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p)
 {
-  G4double tmax = kinEnergy;
-  if(isElectron) tmax *= 0.5;
-  return tmax;
+  particle = p;
+  if(p != theElectron) isElectron = false;
 }
 
Index: trunk/source/processes/electromagnetic/standard/include/G4PAIModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4PAIModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4PAIModel.hh	(revision 1055)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4PAIModel.hh,v 1.22 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -75,12 +77,12 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  virtual void InitialiseMe(const G4ParticleDefinition*) {};
-
-  virtual G4double ComputeDEDX(const G4MaterialCutsCouple*,
+  virtual void InitialiseMe(const G4ParticleDefinition*);
+
+  virtual G4double ComputeDEDXPerVolume(const G4Material*,
 			       const G4ParticleDefinition*,
 			       G4double kineticEnergy,
 			       G4double cutEnergy);
 
-  virtual G4double CrossSection(const G4MaterialCutsCouple*,
+  virtual G4double CrossSectionPerVolume(const G4Material*,
 				const G4ParticleDefinition*,
 				G4double kineticEnergy,
@@ -118,6 +120,4 @@
 
   void SetVerboseLevel(G4int verbose){fVerbose=verbose;};
-
-
 
 protected:
@@ -192,28 +192,4 @@
 };
 
-/////////////////////////////////////////////////////////////////////
-
-inline G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
-                                                      G4double kinEnergy) 
-{
-  G4double tmax = kinEnergy;
-  if(p == fElectron) tmax *= 0.5;
-  else if(p != fPositron) { 
-    G4double mass = p->GetPDGMass();
-    G4double ratio= electron_mass_c2/mass;
-    G4double gamma= kinEnergy/mass + 1.0;
-    tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
-                  (1. + 2.0*gamma*ratio + ratio*ratio);
-  }
-  return tmax;
-}
-
-///////////////////////////////////////////////////////////////
-
-inline  void G4PAIModel::DefineForRegion(const G4Region* r) 
-{
-  fPAIRegionVector.push_back(r);
-}
-
 #endif
 
Index: trunk/source/processes/electromagnetic/standard/include/G4PAIPhotonModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4PAIPhotonModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4PAIPhotonModel.hh	(revision 1055)
@@ -24,4 +24,6 @@
 // ********************************************************************
 //
+// $Id: G4PAIPhotonModel.hh,v 1.12 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -74,16 +76,16 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
   
-  virtual void InitialiseMe(const G4ParticleDefinition*) {};
-
-  virtual G4double ComputeDEDX(const G4MaterialCutsCouple*,
-			       const G4ParticleDefinition*,
-			       G4double kineticEnergy,
-			       G4double cutEnergy);
-
-  virtual G4double CrossSection(const G4MaterialCutsCouple*,
-				const G4ParticleDefinition*,
-				G4double kineticEnergy,
-				G4double cutEnergy,
-				G4double maxEnergy);
+  virtual void InitialiseMe(const G4ParticleDefinition*);
+
+  virtual G4double ComputeDEDXPerVolume(const G4Material*,
+					const G4ParticleDefinition*,
+					G4double kineticEnergy,
+					G4double cutEnergy);
+
+  virtual G4double CrossSectionPerVolume(const G4Material*,
+					 const G4ParticleDefinition*,
+					 G4double kineticEnergy,
+					 G4double cutEnergy,
+					 G4double maxEnergy);
 
   virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
@@ -122,10 +124,8 @@
                              G4double position, G4int iTransfer );
 
-
-
 protected:
 
   G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-                                    G4double kinEnergy);
+			      G4double kinEnergy);
 
 private:
@@ -145,6 +145,4 @@
   G4int                fVerbose; 
   G4PhysicsLogVector*  fProtonEnergyVector ;
-
-
 
   // vectors
@@ -204,28 +202,4 @@
 };
 
-/////////////////////////////////////////////////////////////////////
-
-inline G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
-                                                      G4double kinEnergy) 
-{
-  G4double tmax = kinEnergy;
-  if(p == fElectron) tmax *= 0.5;
-  else if(p != fPositron) { 
-    G4double mass = p->GetPDGMass();
-    G4double ratio= electron_mass_c2/mass;
-    G4double gamma= kinEnergy/mass + 1.0;
-    tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
-                  (1. + 2.0*gamma*ratio + ratio*ratio);
-  }
-  return tmax;
-}
-
-///////////////////////////////////////////////////////////////
-
-inline  void G4PAIPhotonModel::DefineForRegion(const G4Region* r) 
-{
-  fPAIRegionVector.push_back(r);
-}
-
 #endif
 
Index: trunk/source/processes/electromagnetic/standard/include/G4PEEffectModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4PEEffectModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4PEEffectModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4PEEffectModel.hh,v 1.6 2007/05/22 17:34:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4PEEffectModel.hh,v 1.7 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -41,4 +41,6 @@
 //
 // 06.02.2006 : added ComputeMeanFreePath()  (mma)
+// 20.02.2009 : move virtual inline to .cc, substitute
+//              ComputeMeanFreePath() by CrossSectionPerVolume (VI)
 //
 // Class Description:
@@ -70,14 +72,15 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
-                                      G4double kinEnergy,
-                                      G4double Z,
-                                      G4double A,
-                                      G4double, G4double);
+  virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
+					      G4double kinEnergy,
+					      G4double Z,
+					      G4double A,
+					      G4double, G4double);
 				      
-  G4double ComputeMeanFreePath( const G4ParticleDefinition*,
-                                      G4double kinEnergy,
-				const G4Material* material,      
-                                      G4double, G4double);
+  virtual G4double CrossSectionPerVolume(const G4Material*,
+					 const G4ParticleDefinition*,
+					 G4double kineticEnergy,
+					 G4double cutEnergy,
+					 G4double maxEnergy);
 
   virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
@@ -100,43 +103,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
-
-inline G4double G4PEEffectModel::ComputeCrossSectionPerAtom(
-                                       const G4ParticleDefinition*,
-                                             G4double energy,
-                                             G4double Z, G4double,
-                                             G4double, G4double)
-{
- G4double* SandiaCof = G4SandiaTable::GetSandiaCofPerAtom((G4int)Z, energy);
-
- G4double energy2 = energy*energy, energy3 = energy*energy2,
-          energy4 = energy2*energy2;
-
- return SandiaCof[0]/energy  + SandiaCof[1]/energy2 +
-        SandiaCof[2]/energy3 + SandiaCof[3]/energy4;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4double G4PEEffectModel::ComputeMeanFreePath(
-                                       const G4ParticleDefinition*,
-                                             G4double energy,
-				       const G4Material* material,
-                                             G4double, G4double)
-{
- G4double* SandiaCof = material->GetSandiaTable()
-                                ->GetSandiaCofForMaterial(energy);
-				
- G4double energy2 = energy*energy, energy3 = energy*energy2,
-          energy4 = energy2*energy2;
-	  
- G4double cross = SandiaCof[0]/energy  + SandiaCof[1]/energy2 +
-                  SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 
- 
- G4double mfp = DBL_MAX;
- if (cross > 0.) mfp = 1./cross;
- return mfp;
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
Index: trunk/source/processes/electromagnetic/standard/include/G4PhotoElectricEffect.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4PhotoElectricEffect.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4PhotoElectricEffect.hh	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4PhotoElectricEffect.hh,v 1.24 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4PhotoElectricEffect.hh,v 1.25 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //
@@ -91,12 +91,12 @@
 
   // true for Gamma only.
-  G4bool IsApplicable(const G4ParticleDefinition&);
+  virtual G4bool IsApplicable(const G4ParticleDefinition&);
 
   // Print few lines of informations about the process: validity range,
-  void PrintInfo();
+  virtual void PrintInfo();
 
 protected:
 
-  void InitialiseProcess(const G4ParticleDefinition*);
+  virtual void InitialiseProcess(const G4ParticleDefinition*);
 
 private:
@@ -105,13 +105,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline 
-G4bool G4PhotoElectricEffect::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Gamma::Gamma());
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
Index: trunk/source/processes/electromagnetic/standard/include/G4UniversalFluctuation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UniversalFluctuation.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4UniversalFluctuation.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UniversalFluctuation.hh,v 1.6 2008/10/22 16:04:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UniversalFluctuation.hh,v 1.10 2009/03/19 14:15:17 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -70,21 +70,19 @@
   virtual ~G4UniversalFluctuation();
 
-  G4double SampleFluctuations(const G4Material*,
-                          const G4DynamicParticle*,
- 				G4double&,
-                                G4double&,
-                                G4double&);
+  virtual G4double SampleFluctuations(const G4Material*,
+				      const G4DynamicParticle*,
+				      G4double&,
+				      G4double&,
+				      G4double&);
 
-  G4double Dispersion(    const G4Material*,
-                          const G4DynamicParticle*,
- 				G4double&,
-                                G4double&);
+  virtual G4double Dispersion(    const G4Material*,
+				  const G4DynamicParticle*,
+				  G4double&,
+				  G4double&);
 
-  void InitialiseMe(const G4ParticleDefinition*);
+  virtual void InitialiseMe(const G4ParticleDefinition*);
 
   // Initialisation prestep
-  inline void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2);
-
-protected:
+  virtual void SetParticleAndCharge(const G4ParticleDefinition*, G4double q2);
 
 private:
@@ -113,7 +111,4 @@
   G4double e0;
 
-  G4double facwidth;
-  G4double oldloss;
-  G4double samestep;
   G4double e1,e2;
 
@@ -126,17 +121,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline void 
-G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
-					     G4double q2)
-{
-  if(part != particle) {
-    particle       = part;
-    particleMass   = part->GetPDGMass();
-  }
-  chargeSquare = q2;
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
Index: trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UrbanMscModel.hh,v 1.33 2008/03/10 10:39:21 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UrbanMscModel.hh,v 1.35 2009/04/29 13:30:22 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -37,44 +37,10 @@
 // Author:        Laszlo Urban
 //
-// Creation date: 01.03.2001
+// Creation date: 06.03.2008
 //
 // Modifications:
 //
-// 27-03-03 Move model part from G4MultipleScattering (V.Ivanchenko)
-// 27-03-03 Rename (V.Ivanchenko)
-//
-// 05-08-03 angle distribution has been modified (L.Urban)
-// 26-11-03 new data member currentRange (L.Urban)
-// 01-03-04 changes in data members + signature changed in SampleCosineTheta
-// 11-03-04 changes in data members (L.Urban)
-// 23-04-04 changes in data members and in signature of SampleCosineTheta
-//          (L.Urban)
-// 17-08-04 name of data member facxsi changed to factail (L.Urban)
-// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
-// 15-04-05 optimize internal interface - add SampleSecondaries method (V.Ivanchenko)
-// 11-08-05 computation of lateral correlation added (L.Urban)
-// 02-10-05 nuclear size correction computation removed, the correction
-//          included in the (theoretical) tabulated values (L.Urban)
-// 16-02-06 data members b and xsi have been removed (L.Urban)
-// 17-02-06 Save table of transport cross sections not mfp (V.Ivanchenko)
-// 07-03-06 Create G4UrbanMscModel and move there step limit calculation (V.Ivanchenko)
-// 10-05-06 SetMscStepLimitation at initialisation (V.Ivantchenko)
-// 11-05-06 name of data member safety changed to presafety, some new data
-//          members added (frscaling1,frscaling2,tlimitminfix,nstepmax)
-//         (L.Urban)
-// 13-10-06 data member factail removed, data member tkinlimit changed
-//          to lambdalimit,
-//          new data members tgeom,tnow,skin,skindepth,Zeff,geomlimit
-//          G4double GeomLimit(const G4Track& track) changed to
-//              void GeomLimit(const G4Track& track) (L.Urban)
-// 20-10-06 parameter theta0 now computed in the (public)
-//          function ComputeTheta0,
-//          single scattering modified allowing not small
-//          angles as well (L.Urban)
-// 31-01-07 code cleaning (L.Urban)
-// 06-02-07 Move SetMscStepLimitation method into the source (VI)
-// 15-02-07 new data member : smallstep (L.Urban)
-// 10-04-07 remove navigator, smallstep, tnow (V.Ivanchenko)
-//
+// 28-04-09  move G4UrbanMscModel2 from the g49.2 to G4UrbanMscModel. 
+//           now it is frozen (V.Ivanchenk0)
 //
 // Class Description:
@@ -111,5 +77,5 @@
 
   void Initialise(const G4ParticleDefinition*, const G4DataVector&);
-				          
+
   G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
 				      G4double KineticEnergy,
@@ -119,10 +85,4 @@
 				      G4double emax=DBL_MAX);
 
-  void SampleSecondaries(std::vector<G4DynamicParticle*>*, 
-			 const G4MaterialCutsCouple*,
-			 const G4DynamicParticle*,
-			 G4double,
-			 G4double);
-
   void SampleScattering(const G4DynamicParticle*,
 			G4double safety);
@@ -141,4 +101,6 @@
 private:
 
+  G4double SimpleScattering(G4double xmeanth, G4double x2meanth);
+
   G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
 
@@ -147,9 +109,9 @@
   G4double LatCorrelation();
 
-  void GeomLimit(const G4Track& track);
-
   inline G4double GetLambda(G4double kinEnergy);
 
   inline void SetParticle(const G4ParticleDefinition*);
+
+  inline void UpdateCache();
 
   //  hide assignment operator
@@ -160,5 +122,4 @@
   G4ParticleChangeForMSC*     fParticleChange;
 
-  G4SafetyHelper*             safetyHelper;
   G4PhysicsTable*             theLambdaTable;
   const G4MaterialCutsCouple* couple;
@@ -166,7 +127,6 @@
 
   G4double mass;
-  G4double charge;
-
-  G4double masslimite,masslimitmu;
+  G4double charge,ChargeSquare;
+  G4double masslimite,lambdalimit,fr;
 
   G4double taubig;
@@ -174,10 +134,9 @@
   G4double taulim;
   G4double currentTau;
-  G4double frscaling1,frscaling2;
   G4double tlimit;
   G4double tlimitmin;
   G4double tlimitminfix;
-
-  G4double nstepmax;
+  G4double tgeom;
+
   G4double geombig;
   G4double geommin;
@@ -198,16 +157,25 @@
   G4double currentKinEnergy;
   G4double currentRange; 
+  G4double rangeinit;
   G4double currentRadLength;
 
-  G4double Zeff;
+  G4double theta0max,rellossmax;
+  G4double third;
 
   G4int    currentMaterialIndex;
+
+  G4double y;
+  G4double Zold;
+  G4double Zeff,Z2,Z23,lnZ;
+  G4double coeffth1,coeffth2;
+  G4double coeffc1,coeffc2;
+  G4double scr1ini,scr2ini,scr1,scr2;
 
   G4bool   isInitialized;
   G4bool   inside;
   G4bool   insideskin;
-
 };
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -236,4 +204,5 @@
     mass = p->GetPDGMass();
     charge = p->GetPDGCharge()/eplus;
+    ChargeSquare = charge*charge;
   }
 }
@@ -241,4 +210,20 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+inline
+void G4UrbanMscModel::UpdateCache()                                   
+{
+    lnZ = std::log(Zeff);
+    coeffth1 = 0.885+lnZ*(0.104-0.0170*lnZ);
+    coeffth2 = 0.028+lnZ*(0.012-0.00125*lnZ);
+    coeffc1  = 2.134-lnZ*(0.1045-0.00602*lnZ);
+    coeffc2  = 0.001126-lnZ*(0.0001089+0.0000247*lnZ);
+    Z2 = Zeff*Zeff;
+    Z23 = std::exp(2.*lnZ/3.);
+    scr1     = scr1ini*Z23;
+    scr2     = scr2ini*Z2*ChargeSquare;
+  //  lastMaterial = couple->GetMaterial();
+    Zold = Zeff;
+}
+
 #endif
 
Index: trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel2.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel2.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel2.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UrbanMscModel2.hh,v 1.11 2008/12/18 13:01:34 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UrbanMscModel2.hh,v 1.17 2009/05/15 09:26:42 urban Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -40,5 +40,5 @@
 //
 // Modifications:
-//
+// 23.04.2009 L.Urban updated parameterization in UpdateCache method
 //
 //
@@ -84,10 +84,4 @@
 				      G4double emax=DBL_MAX);
 
-  void SampleSecondaries(std::vector<G4DynamicParticle*>*, 
-			 const G4MaterialCutsCouple*,
-			 const G4DynamicParticle*,
-			 G4double,
-			 G4double);
-
   void SampleScattering(const G4DynamicParticle*,
 			G4double safety);
@@ -114,6 +108,4 @@
   G4double LatCorrelation();
 
-  void GeomLimit(const G4Track& track);
-
   inline G4double GetLambda(G4double kinEnergy);
 
@@ -129,9 +121,7 @@
   G4ParticleChangeForMSC*     fParticleChange;
 
-  G4SafetyHelper*             safetyHelper;
   G4PhysicsTable*             theLambdaTable;
   const G4MaterialCutsCouple* couple;
   G4LossTableManager*         theManager;
-
 
   G4double mass;
@@ -223,13 +213,17 @@
 {
     lnZ = std::log(Zeff);
-    coeffth1 = 0.885+lnZ*(0.104-0.0170*lnZ);
-    coeffth2 = 0.028+lnZ*(0.012-0.00125*lnZ);
-    coeffc1  = 2.134-lnZ*(0.1045-0.00602*lnZ);
-    coeffc2  = 0.001126-lnZ*(0.0001089+0.0000247*lnZ);
+    //new correction in theta0 formula
+    coeffth1 = (1.-8.7780e-2/Zeff)*(0.87+0.03*lnZ);                   
+    coeffth2 = (4.0780e-2+1.7315e-4*Zeff)*(0.87+0.03*lnZ);              
+    // tail parameters
+    G4double lnZ1 = std::log(Zeff+1.);
+    coeffc1  = 2.943-0.197*lnZ1;                  
+    coeffc2  = 0.0987-0.0143*lnZ1;                              
+    // for single scattering
     Z2 = Zeff*Zeff;
     Z23 = std::exp(2.*lnZ/3.);
     scr1     = scr1ini*Z23;
     scr2     = scr2ini*Z2*ChargeSquare;
-  //  lastMaterial = couple->GetMaterial();
+                                              
     Zold = Zeff;
 }
Index: trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel90.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel90.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel90.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UrbanMscModel90.hh,v 1.4 2008/10/29 14:15:30 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UrbanMscModel90.hh,v 1.5 2009/04/10 16:34:56 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -110,6 +110,4 @@
 
   G4double LatCorrelation();
-
-  void GeomLimit(const G4Track& track);
 
   inline G4double GetLambda(G4double kinEnergy);
Index: trunk/source/processes/electromagnetic/standard/include/G4WentzelVIModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4WentzelVIModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4WentzelVIModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WentzelVIModel.hh,v 1.7 2008/08/04 08:49:09 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4WentzelVIModel.hh,v 1.17 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -66,5 +66,4 @@
 class G4LossTableManager;
 class G4ParticleChangeForMSC;
-class G4SafetyHelper;
 class G4ParticleDefinition;
 
@@ -80,32 +79,26 @@
   virtual ~G4WentzelVIModel();
 
-  void Initialise(const G4ParticleDefinition*, const G4DataVector&);
-
-  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
-				      G4double KineticEnergy,
-				      G4double AtomicNumber,
-				      G4double AtomicWeight=0., 
-				      G4double cut = DBL_MAX,
-				      G4double emax= DBL_MAX);
-
-  void SampleScattering(const G4DynamicParticle*, G4double safety);
-
-  void SampleSecondaries(std::vector<G4DynamicParticle*>*, 
-			 const G4MaterialCutsCouple*,
-			 const G4DynamicParticle*,
-			 G4double,
-			 G4double);
-
-  G4double ComputeTruePathLengthLimit(const G4Track& track,
-				      G4PhysicsTable* theLambdaTable,
-				      G4double currentMinimalStep);
-
-  G4double ComputeGeomPathLength(G4double truePathLength);
-
-  G4double ComputeTrueStepLength(G4double geomStepLength);
+  virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
+
+  virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
+					      G4double KineticEnergy,
+					      G4double AtomicNumber,
+					      G4double AtomicWeight=0., 
+					      G4double cut = DBL_MAX,
+					      G4double emax= DBL_MAX);
+
+  virtual void SampleScattering(const G4DynamicParticle*, G4double safety);
+
+  virtual G4double ComputeTruePathLengthLimit(const G4Track& track,
+					      G4PhysicsTable* theLambdaTable,
+					      G4double currentMinimalStep);
+
+  virtual G4double ComputeGeomPathLength(G4double truePathLength);
+
+  virtual G4double ComputeTrueStepLength(G4double geomStepLength);
 
 private:
 
-  G4double ComputeTransportXSectionPerVolume();
+  G4double ComputeTransportXSectionPerAtom();
 
   G4double ComputeXSectionPerVolume();
@@ -133,5 +126,4 @@
   G4ParticleChangeForMSC*   fParticleChange;
 
-  G4SafetyHelper*           safetyHelper;
   G4PhysicsTable*           theLambdaTable;
   G4PhysicsTable*           theLambda2Table;
@@ -173,5 +165,4 @@
   // single scattering parameters
   G4double coeff;
-  G4double constn;
   G4double cosThetaMin;
   G4double cosThetaMax;
@@ -182,5 +173,4 @@
   G4double q2Limit;
   G4double alpha2;
-  G4double a0;
 
   // projectile
@@ -193,4 +183,5 @@
   G4double mom2;
   G4double invbeta2;
+  G4double kinFactor;
   G4double etag;
   G4double lowEnergyLimit;
@@ -198,7 +189,11 @@
   // target
   G4double targetZ;
+  G4double targetMass;
   G4double screenZ;
   G4double formfactA;
-  G4double FF[100];
+  G4int    iz;
+
+  static G4double ScreenRSquare[100];
+  static G4double FormFactor[100];
 
   // flags
@@ -240,6 +235,5 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline 
-void G4WentzelVIModel::SetupParticle(const G4ParticleDefinition* p)
+inline void G4WentzelVIModel::SetupParticle(const G4ParticleDefinition* p)
 {
   // Initialise mass and charge
@@ -251,5 +245,4 @@
     chargeSquare = q*q;
     tkin = 0.0;
-    lowEnergyLimit = keV*mass/electron_mass_c2;
   }
 }
@@ -264,5 +257,5 @@
     invbeta2 = 1.0 +  mass*mass/mom2;
     cosTetMaxNuc = cosThetaMax;
-    if(ekin <= 10.*cut && mass < MeV) {
+    if(mass < MeV && ekin <= 10.*cut) {
       cosTetMaxNuc = ekin*(cosThetaMax + 1.0)/(10.*cut) - 1.0;
     }
@@ -278,26 +271,18 @@
     etag    = e; 
     targetZ = Z;
-    G4int iz= G4int(Z);
+    iz = G4int(Z);
     if(iz > 99) iz = 99;
-    G4double x = fNistManager->GetZ13(iz);
-    screenZ = a0*x*x/mom2;
-    if(iz > 1) screenZ *=(1.13 + 3.76*invbeta2*Z*Z*chargeSquare*alpha2);
-    //    screenZ = a0*x*x*(1.13 + 3.76*Z*Z*chargeSquare*alpha2)/mom2;
-    // A.V. Butkevich et al., NIM A 488 (2002) 282
-    formfactA = FF[iz];
-    if(formfactA == 0.0) {
-      x = fNistManager->GetA27(iz); 
-      formfactA = constn*x*x;
-      FF[iz] = formfactA;
+    targetMass = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
+    G4double m12  = mass*mass;
+    G4double x    = 1.0 + mass/targetMass;
+    kinFactor  = coeff*targetZ*chargeSquare*(1.0 +  m12/mom2)/mom2;
+    screenZ = ScreenRSquare[iz]/mom2;
+    if(iz > 1) {
+      screenZ *=(1.13 + 3.76*Z*Z*alpha2);
+      kinFactor /= (x*x);
     }
-    formfactA *= mom2;
+    //if(iz > 1) screenZ *=(1.13 + std::min(0.5,3.76*Z*Z*invbeta2*alpha2));
+    formfactA = FormFactor[iz]*mom2;
     cosTetMaxNuc2 = cosTetMaxNuc;
-    /*
-    G4double ee = 10.*eV*Z;
-    if(1 == iz) ee *= 2.0;
-    G4double z = std::min(cosTetMaxElec, 1.0 - std::max(ecut,ee)*amu_c2
-			  *fNistManager->GetAtomicMassAmu(iz)/mom2);
-    cosTetMaxElec2 = std::max(cosTetMaxNuc2, z);
-    */
     cosTetMaxElec2 = cosTetMaxElec;
   } 
Index: trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlung.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlung.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlung.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlung.hh,v 1.36 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlung.hh,v 1.37 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -89,5 +89,5 @@
   virtual ~G4eBremsstrahlung();
 
-  G4bool IsApplicable(const G4ParticleDefinition& p);
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
 
   // Print out of the class parameters
@@ -111,12 +111,4 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4eBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 #endif
Index: trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlungModel.hh,v 1.25 2008/11/13 19:28:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlungModel.hh,v 1.26 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -75,6 +75,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*, 
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*, 
+				const G4MaterialCutsCouple*);
 
   virtual G4double ComputeDEDXPerVolume(const G4Material*,
@@ -102,7 +102,4 @@
 
 protected:
-
-  inline G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-				     G4double kineticEnergy);
 
   const G4Element* SelectRandomAtom(const G4MaterialCutsCouple* couple);
@@ -191,13 +188,3 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline 
-G4double G4eBremsstrahlungModel::MaxSecondaryEnergy(
-                                 const G4ParticleDefinition*,
-    				       G4double kineticEnergy)
-{
-  return kineticEnergy;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 #endif
Index: trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungRelModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungRelModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eBremsstrahlungRelModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlungRelModel.hh,v 1.10 2008/11/14 09:25:19 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlungRelModel.hh,v 1.11 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -71,6 +71,6 @@
   virtual void Initialise(const G4ParticleDefinition*, const G4DataVector&);
 
-  G4double MinEnergyCut(const G4ParticleDefinition*, 
-			const G4MaterialCutsCouple*);
+  virtual G4double MinEnergyCut(const G4ParticleDefinition*, 
+				const G4MaterialCutsCouple*);
 
   virtual G4double ComputeDEDXPerVolume(const G4Material*,
@@ -97,9 +97,4 @@
   inline G4double LPMconstant() const;
 
-protected:
-
-  inline G4double MaxSecondaryEnergy(const G4ParticleDefinition*,
-				     G4double kineticEnergy);
-
 private:
 
@@ -207,13 +202,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-inline G4double 
-G4eBremsstrahlungRelModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
-					      G4double kineticEnergy)
-{
-  return kineticEnergy;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 
 inline G4double G4eBremsstrahlungRelModel::Phi1(G4double gg, G4double)
Index: trunk/source/processes/electromagnetic/standard/include/G4eCoulombScatteringModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eCoulombScatteringModel.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eCoulombScatteringModel.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eCoulombScatteringModel.hh,v 1.36 2008/08/04 08:49:09 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eCoulombScatteringModel.hh,v 1.43 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -65,6 +65,6 @@
 #include "G4VEmModel.hh"
 #include "G4PhysicsTable.hh"
+#include "globals.hh"
 #include "G4NistManager.hh"
-#include "globals.hh"
 
 class G4ParticleChangeForGamma;
@@ -137,5 +137,4 @@
 
   G4double                  coeff;
-  G4double                  constn;
   G4double                  cosThetaMin;
   G4double                  cosThetaMax;
@@ -160,4 +159,5 @@
   G4double                  mom2;
   G4double                  invbeta2;
+  G4double                  kinFactor;
   G4double                  etag;
   G4double                  lowEnergyLimit;
@@ -165,14 +165,17 @@
   // target
   G4double                  targetZ;
+  G4double                  targetMass;
   G4double                  screenZ;
   G4double                  formfactA;
   G4int                     idxelm;
+  G4int                     iz;
 
 private:
 
-  G4double                  a0;
   G4double                  alpha2;
   G4double                  faclim;
-  G4double                  FF[100];
+
+  static G4double ScreenRSquare[100];
+  static G4double FormFactor[100];
 
   G4bool                    isInitialised;             
@@ -204,5 +207,4 @@
     chargeSquare = q*q;
     tkin = 0.0;
-    lowEnergyLimit = keV*mass/electron_mass_c2;
   }
 }
@@ -219,5 +221,5 @@
     cosTetMinNuc = cosThetaMin;
     cosTetMaxNuc = cosThetaMax;
-    if(ekin <= 10.*cut && mass < MeV && cosThetaMin < 1.0) {
+    if(mass < MeV && cosThetaMin < 1.0 && ekin <= 10.*cut) {
       cosTetMinNuc = ekin*(cosThetaMin + 1.0)/(10.*cut) - 1.0;
     }
@@ -233,29 +235,22 @@
     etag    = e; 
     targetZ = Z;
-    G4int iz= G4int(Z);
+    iz= G4int(Z);
     if(iz > 99) iz = 99;
-    G4double x = fNistManager->GetZ13(iz);
-    screenZ = a0*x*x/mom2;
-    if(iz > 1) screenZ *=(1.13 + 3.76*invbeta2*Z*Z*chargeSquare*alpha2);
-    //screenZ = a0*x*x*(1.13 + 3.76*Z*Z*chargeSquare*alpha2)/mom2;
-    // A.V. Butkevich et al., NIM A 488 (2002) 282
-    formfactA = FF[iz];
-    if(formfactA == 0.0) {
-      x = fNistManager->GetA27(iz); 
-      formfactA = constn*x*x;
-      FF[iz] = formfactA;
+    targetMass   = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
+    G4double m12 = mass*mass;
+    G4double x   = 1.0 + mass/targetMass;
+    kinFactor    = (1.0 +  m12/mom2)/mom2;
+
+    screenZ = ScreenRSquare[iz]/mom2;
+    if(iz > 1) {
+      screenZ *=(1.13 + 3.76*Z*Z*alpha2);
+      kinFactor /= (x*x);
     }
-    formfactA *= mom2;
+    // if(iz > 1) screenZ *=(1.13 + std::min(0.5,3.76*Z*Z*invbeta2*alpha2));
+    formfactA = FormFactor[iz]*mom2;
     cosTetMaxNuc2 = cosTetMaxNuc;
-    if(particle == theProton && 1 == iz && cosTetMaxNuc2 < 0.0) {
+    if(1 == iz && particle == theProton && cosTetMaxNuc2 < 0.0) {
       cosTetMaxNuc2 = 0.0;
     }
-    /*
-    G4double ee = 10.*eV*Z;
-    if(1 == iz) ee *= 2.0;
-    G4double z = std::min(cosTetMaxElec, 1.0 - std::max(ecut,ee)*amu_c2
-			  *fNistManager->GetAtomicMassAmu(iz)/mom2);
-    cosTetMaxElec2 = std::max(cosTetMaxNuc2, z);
-    */
     cosTetMaxElec2 = cosTetMaxElec;
   } 
Index: trunk/source/processes/electromagnetic/standard/include/G4eIonisation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eIonisation.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eIonisation.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eIonisation.hh,v 1.35 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eIonisation.hh,v 1.36 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -88,5 +88,5 @@
   virtual ~G4eIonisation();
 
-  G4bool IsApplicable(const G4ParticleDefinition& p);
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
 
   // Print out of the class parameters
@@ -98,6 +98,6 @@
 					   const G4ParticleDefinition*);
 
-  G4double MinPrimaryEnergy(const G4ParticleDefinition*,
-			    const G4Material*, G4double cut);
+  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition*,
+				    const G4Material*, G4double cut);
 
 private:
@@ -114,23 +114,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4double G4eIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
-						const G4Material*,
-						G4double cut)
-{
-  G4double x = cut;
-  if(isElectron) x += cut;
-  return x;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4eIonisation::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
Index: trunk/source/processes/electromagnetic/standard/include/G4eplusAnnihilation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4eplusAnnihilation.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4eplusAnnihilation.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eplusAnnihilation.hh,v 1.23 2007/05/23 08:47:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eplusAnnihilation.hh,v 1.24 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -73,5 +73,5 @@
   virtual ~G4eplusAnnihilation();
 
-  G4bool IsApplicable(const G4ParticleDefinition& p);
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
 
   virtual G4VParticleChange* AtRestDoIt(
@@ -79,5 +79,5 @@
                              const G4Step& stepData);
 
-  G4double AtRestGetPhysicalInteractionLength(
+  virtual G4double AtRestGetPhysicalInteractionLength(
                              const G4Track& track,
                              G4ForceCondition* condition
@@ -96,22 +96,4 @@
 };
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4eplusAnnihilation::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (&p == G4Positron::Positron());
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline 
-G4double G4eplusAnnihilation::AtRestGetPhysicalInteractionLength(
-                              const G4Track&, G4ForceCondition* condition)
-{
-  *condition = NotForced;
-  return 0.0;
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
Index: trunk/source/processes/electromagnetic/standard/include/G4hIonisation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4hIonisation.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4hIonisation.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4hIonisation.hh,v 1.41 2008/09/14 17:11:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4hIonisation.hh,v 1.42 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -93,8 +93,8 @@
   virtual ~G4hIonisation();
 
-  G4bool IsApplicable(const G4ParticleDefinition& p);
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
 
-  G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
-			    const G4Material*, G4double cut);
+  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
+				    const G4Material*, G4double cut);
 
   // Print out of the class parameters
@@ -123,24 +123,4 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4hIonisation::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV &&
-	 !p.IsShortLived());
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4double G4hIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
-						const G4Material*,
-						G4double cut)
-{
-  G4double x = 0.5*cut/electron_mass_c2;
-  G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
-  return mass*(g - 1.0);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
 inline void G4hIonisation::ActivateNuclearStopping(G4bool val)
Index: trunk/source/processes/electromagnetic/standard/include/G4ionIonisation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4ionIonisation.hh	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/include/G4ionIonisation.hh	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4ionIonisation.hh,v 1.56 2008/09/14 17:11:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ionIonisation.hh,v 1.57 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -85,5 +85,5 @@
   virtual ~G4ionIonisation();
 
-  inline G4bool IsApplicable(const G4ParticleDefinition& p);
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
 
   // Print out of the class parameters
@@ -102,5 +102,5 @@
 					   const G4ParticleDefinition*);
 
-  inline G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
+  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
 				   const G4Material*, G4double cut);
 
@@ -127,21 +127,4 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4bool G4ionIonisation::IsApplicable(const G4ParticleDefinition& p)
-{
-  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived() &&
-          p.GetParticleType() == "nucleus");
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-inline G4double G4ionIonisation::MinPrimaryEnergy(
-          const G4ParticleDefinition* p, const G4Material*, G4double cut)
-{
-  return 
-    p->GetPDGMass()*(std::sqrt(1. + 0.5*cut/CLHEP::electron_mass_c2) - 1.0);
-}
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
Index: trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4BetheBlochModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BetheBlochModel.cc,v 1.24 2008/10/22 16:00:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BetheBlochModel.cc,v 1.31 2009/04/24 14:24:00 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -65,5 +65,4 @@
 #include "G4EmCorrections.hh"
 #include "G4ParticleChangeForLoss.hh"
-#include "G4NistManager.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -83,5 +82,8 @@
 {
   fParticleChange = 0;
-  if(p) SetParticle(p);
+  if(p) {
+    SetGenericIon(p);
+    SetParticle(p);
+  }
   theElectron = G4Electron::Electron();
   corr = G4LossTableManager::Instance()->EmCorrections();  
@@ -108,53 +110,18 @@
                                    const G4DataVector&)
 {
-  if (!particle) SetParticle(p);
+  SetGenericIon(p);
+  SetParticle(p);
+
+  //G4cout << "G4BetheBlochModel::Initialise for " << p->GetParticleName()
+  //	 << "  isIon= " << isIon 
+  //	 << G4endl;
 
   corrFactor = chargeSquare;
+  // always false before the run
+  SetDeexcitationFlag(false);
 
   if(!isInitialised) {
     isInitialised = true;
-
-    if(!fParticleChange) {
-      if (pParticleChange) {
-	fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
-	  (pParticleChange);
-      } else { 
-	fParticleChange = new G4ParticleChangeForLoss();
-      }
-    }
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4BetheBlochModel::SetParticle(const G4ParticleDefinition* p)
-{
-  if(particle != p) {
-    particle = p;
-    G4String pname = particle->GetParticleName();
-    if (particle->GetParticleType() == "nucleus" &&
-	pname != "deuteron" && pname != "triton") {
-      isIon = true;
-    }
-    
-    mass = particle->GetPDGMass();
-    spin = particle->GetPDGSpin();
-    G4double q = particle->GetPDGCharge()/eplus;
-    chargeSquare = q*q;
-    ratio = electron_mass_c2/mass;
-    G4double magmom = particle->GetPDGMagneticMoment()
-      *mass/(0.5*eplus*hbar_Planck*c_squared);
-    magMoment2 = magmom*magmom - 1.0;
-    formfact = 0.0;
-    if(particle->GetLeptonNumber() == 0) {
-      G4double x = 0.8426*GeV;
-      if(spin == 0.0 && mass < GeV) {x = 0.736*GeV;}
-      else if(mass > GeV) {
-        x /= nist->GetZ13(mass/proton_mass_c2);
-	//	tlimit = 51.2*GeV*A13[iz]*A13[iz];
-      }
-      formfact = 2.0*electron_mass_c2/(x*x);
-      tlimit   = 2.0/formfact;
-    }
+    fParticleChange = GetParticleChangeForLoss();
   }
 }
@@ -415,4 +382,5 @@
   G4double cost = deltaKinEnergy * (totEnergy + electron_mass_c2) /
                                    (deltaMomentum * totMomentum);
+  /*
   if(cost > 1.0) {
     G4cout << "### G4BetheBlochModel WARNING: cost= " 
@@ -429,4 +397,5 @@
     cost = 1.0;
   }
+  */
   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
 
@@ -440,5 +409,5 @@
   // create G4DynamicParticle object for delta ray
   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
-                                                 deltaDirection,deltaKinEnergy);
+						   deltaDirection,deltaKinEnergy);
 
   vdp->push_back(delta);
@@ -453,3 +422,16 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4BetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
+					       G4double kinEnergy) 
+{
+  // here particle type is checked for any method
+  SetParticle(pd);
+  G4double tau  = kinEnergy/mass;
+  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
+                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
+  return std::min(tmax,tlimit);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4BetheHeitlerModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BetheHeitlerModel.cc,v 1.12 2008/10/15 15:54:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BetheHeitlerModel.cc,v 1.13 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -94,11 +94,5 @@
 				     const G4DataVector&)
 {
-  if(!fParticleChange) {
-    if(pParticleChange) {
-      fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-    } else {
-      fParticleChange = new G4ParticleChangeForGamma();
-    }
-  }
+  if(!fParticleChange) fParticleChange = GetParticleChangeForGamma();
 
   if(theCrossSectionTable) {
Index: trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4BohrFluctuations.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BohrFluctuations.cc,v 1.6 2007/09/27 14:02:41 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BohrFluctuations.cc,v 1.7 2009/02/19 19:17:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -137,3 +137,22 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4double G4BohrFluctuations::Dispersion(const G4Material* material,
+					const G4DynamicParticle* dp,
+					G4double& tmax,
+					G4double& length)
+{
+  if(!particle) InitialiseMe(dp->GetDefinition());
 
+  G4double electronDensity = material->GetElectronDensity();
+  kineticEnergy = dp->GetKineticEnergy();
+  G4double etot = kineticEnergy + particleMass;
+  beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
+  G4double siga  = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * length
+                 * electronDensity * chargeSquare;
+
+  return siga;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
Index: trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4BraggIonModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggIonModel.cc,v 1.22 2008/10/22 16:00:57 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BraggIonModel.cc,v 1.24 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -114,4 +114,7 @@
   corrFactor = chargeSquare;
 
+  // always false before the run
+  SetDeexcitationFlag(false);
+
   if(!isInitialised) {
     isInitialised = true;
@@ -123,12 +126,5 @@
     corr = G4LossTableManager::Instance()->EmCorrections();
 
-    if(!fParticleChange) {
-      if(pParticleChange) {
-	fParticleChange = 
-	  reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-      } else {
-	fParticleChange = new G4ParticleChangeForLoss();
-      }
-    }
+    fParticleChange = GetParticleChangeForLoss();
   }
 }
@@ -352,4 +348,16 @@
   fParticleChange->SetProposedKineticEnergy(kineticEnergy);
   fParticleChange->SetProposedMomentumDirection(finalP);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BraggIonModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
+					     G4double kinEnergy)
+{
+  if(pd != particle) SetParticle(pd);
+  G4double tau  = kinEnergy/mass;
+  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
+                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
+  return tmax;
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4BraggModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4BraggModel.cc,v 1.20 2008/10/22 16:01:46 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4BraggModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -112,4 +112,7 @@
   if(p != particle) SetParticle(p);
 
+  // always false before the run
+  SetDeexcitationFlag(false);
+
   if(!isInitialised) {
     isInitialised = true;
@@ -121,10 +124,5 @@
     corr = G4LossTableManager::Instance()->EmCorrections();
 
-    if(pParticleChange) {
-      fParticleChange = 
-	reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-    } else {
-      fParticleChange = new G4ParticleChangeForLoss();
-    }
+    fParticleChange = GetParticleChangeForLoss();
   }
 }
@@ -337,4 +335,16 @@
 
   vdp->push_back(delta);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4BraggModel::MaxSecondaryEnergy(const G4ParticleDefinition* pd,
+					  G4double kinEnergy)
+{
+  if(pd != particle) SetParticle(pd);
+  G4double tau  = kinEnergy/mass;
+  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
+                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
+  return tmax;
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4ComptonScattering.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4ComptonScattering.cc,v 1.30 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ComptonScattering.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // 
@@ -84,4 +84,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4bool G4ComptonScattering::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Gamma::Gamma());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4ComptonScattering::InitialiseProcess(const G4ParticleDefinition*)
 {
Index: trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4CoulombScattering.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScattering.cc,v 1.19 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4CoulombScattering.cc,v 1.20 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -79,4 +79,11 @@
 G4CoulombScattering::~G4CoulombScattering()
 {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4CoulombScattering::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived());
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4CoulombScatteringModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4CoulombScatteringModel.cc,v 1.37 2008/07/31 13:11:34 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4CoulombScatteringModel.cc,v 1.39 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -93,5 +93,5 @@
 
   // CM system
-  G4int iz      = G4int(Z);
+  iz            = G4int(Z);
   G4double m2   = fNistManager->GetAtomicMassAmu(iz)*amu_c2;
   G4double etot = tkin + mass;
@@ -137,5 +137,5 @@
 
   G4double Z  = currentElement->GetZ();
-  G4int iz    = G4int(Z);
+  iz          = G4int(Z);
   G4int ia    = SelectIsotopeNumber(currentElement);
   G4double m2 = theParticleTable->GetIonTable()->GetNucleusMass(iz, ia);
@@ -182,9 +182,6 @@
   // recoil
   G4double erec = kinEnergy - ekin;
-  G4double th = 
-    std::min(recoilThreshold,
-	     Z*currentElement->GetIonisation()->GetMeanExcitationEnergy());
 
-  if(erec > th) {
+  if(erec > recoilThreshold) {
     G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
     G4double plab = sqrt(ekin*(ekin + 2.0*mass));
Index: trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4GammaConversion.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4GammaConversion.cc,v 1.30 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4GammaConversion.cc,v 1.31 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // 
@@ -90,4 +90,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4bool G4GammaConversion::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Gamma::Gamma());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4GammaConversion::InitialiseProcess(const G4ParticleDefinition*)
 {
Index: trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4IonFluctuations.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4IonFluctuations.cc,v 1.24 2008/10/22 16:25:21 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4IonFluctuations.cc,v 1.26 2009/03/31 13:24:40 toshito Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -190,14 +190,16 @@
 
   // heavy ion correction
-  G4double f1 = 1.065e-4*chargeSquare;
-  if(beta2 > theBohrBeta2)  f1/= beta2;
-  else                      f1/= theBohrBeta2;
-  if(f1 > 2.5) f1 = 2.5;
-  fac *= (1.0 + f1);
+//  G4double f1 = 1.065e-4*chargeSquare;
+//  if(beta2 > theBohrBeta2)  f1/= beta2;
+//  else                      f1/= theBohrBeta2;
+//  if(f1 > 2.5) f1 = 2.5;
+//  fac *= (1.0 + f1);
 
   // taking into account the cut
-  if(fac > 1.0) {
-    siga *= (1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2)));
-  }
+  G4double fac_cut = 1.0 + (fac - 1.0)*2.0*electron_mass_c2*beta2/(tmax*(1.0 - beta2));
+  if(fac_cut > 0.01 && fac > 0.01) {
+    siga *= fac_cut;
+  }
+
   //G4cout << "siga(keV)= " << sqrt(siga)/keV << " fac= " << fac 
   //	 << "  f1= " << f1 << G4endl;
@@ -420,2 +422,17 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4IonFluctuations::SetParticleAndCharge(const G4ParticleDefinition* part,
+					     G4double q2)
+{
+  if(part != particle) {
+    particle       = part;
+    particleMass   = part->GetPDGMass();
+    charge         = part->GetPDGCharge()/eplus;
+    chargeSquare   = charge*charge;
+  }
+  effChargeSquare  = q2;
+  uniFluct.SetParticleAndCharge(part, q2);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4KleinNishinaCompton.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4KleinNishinaCompton.cc,v 1.9 2007/05/22 17:34:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4KleinNishinaCompton.cc,v 1.10 2009/05/15 17:12:33 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -79,8 +79,5 @@
                                        const G4DataVector&)
 {
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForGamma();
+  fParticleChange = GetParticleChangeForGamma();
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4MollerBhabhaModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MollerBhabhaModel.cc,v 1.30 2007/05/22 17:34:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4MollerBhabhaModel.cc,v 1.34 2009/04/24 17:15:46 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -74,8 +74,9 @@
                                          const G4String& nam)
   : G4VEmModel(nam),
-  particle(0),
-  isElectron(true),
-  twoln10(2.0*log(10.0)),
-  lowLimit(0.2*keV)
+    particle(0),
+    isElectron(true),
+    twoln10(2.0*log(10.0)),
+    lowLimit(0.2*keV),
+    isInitialised(false)
 {
   theElectron = G4Electron::Electron();
@@ -87,12 +88,4 @@
 G4MollerBhabhaModel::~G4MollerBhabhaModel()
 {}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MollerBhabhaModel::SetParticle(const G4ParticleDefinition* p)
-{
-  particle = p;
-  if(p != theElectron) isElectron = false;
-}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -108,13 +101,24 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
+G4double G4MollerBhabhaModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
+						 G4double kinEnergy) 
+{
+  G4double tmax = kinEnergy;
+  if(isElectron) tmax *= 0.5;
+  return tmax;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 void G4MollerBhabhaModel::Initialise(const G4ParticleDefinition* p,
                                      const G4DataVector&)
 {
   if(!particle) SetParticle(p);
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
-                                                     (pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+  SetDeexcitationFlag(false);
+
+  if(isInitialised) return;
+
+  isInitialised = true;
+  fParticleChange = GetParticleChangeForLoss();
 }
 
@@ -289,8 +293,10 @@
 					    G4double maxEnergy)
 {
-  G4double tmax = std::min(maxEnergy, MaxSecondaryKinEnergy(dp));
+  G4double kineticEnergy = dp->GetKineticEnergy();
+  G4double tmax = kineticEnergy;
+  if(isElectron) tmax *= 0.5;
+  if(maxEnergy < tmax) tmax = maxEnergy;
   if(tmin >= tmax) return;
 
-  G4double kineticEnergy = dp->GetKineticEnergy();
   G4double energy = kineticEnergy + electron_mass_c2;
   G4double totalMomentum = sqrt(kineticEnergy*(energy + electron_mass_c2));
@@ -339,23 +345,11 @@
     G4double b3  = b4 + y122;
 
-    y     = xmax*xmax;
-    grej  = -xmin*b1;
-    grej += y*b2;
-    grej -= xmin*xmin*xmin*b3;
-    grej += y*y*b4;
-    grej *= beta2;
-    grej += 1.0;
+    y    = xmax*xmax;
+    grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + y*b2 - xmin*b1)*beta2; 
     do {
-      q  = G4UniformRand();
-      x  = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
-      z  = -x*b1;
-      y  = x*x;
-      z += y*b2;
-      y *= x;
-      z -= y*b3;
-      y *= x;
-      z += y*b4;
-      z *= beta2;
-      z += 1.0;
+      q = G4UniformRand();
+      x = xmin*xmax/(xmin*(1.0 - q) + xmax*q);
+      y = x*x;
+      z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2; 
       /*
       if(z > grej) {
@@ -396,5 +390,5 @@
   // create G4DynamicParticle object for delta ray
   G4DynamicParticle* delta = new G4DynamicParticle(theElectron,
-						 deltaDirection,deltaKinEnergy);
+						   deltaDirection,deltaKinEnergy);
   vdp->push_back(delta);
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4MscModel71.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MscModel71.cc,v 1.6 2008/03/13 17:20:07 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4MscModel71.cc,v 1.7 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -140,5 +140,5 @@
 
   if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
+    fParticleChange = static_cast<G4ParticleChangeForMSC*>(pParticleChange);
   else
     fParticleChange = new G4ParticleChangeForMSC();
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIModel.cc	(revision 1055)
@@ -24,4 +24,10 @@
 // ********************************************************************
 //
+// $Id: G4PAIModel.cc,v 1.47 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class
 // File name:     G4PAIModel.cc
 //
@@ -153,8 +159,5 @@
 						 fTotBin);
 
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+  fParticleChange = GetParticleChangeForLoss();
 
   // Prepare initialization
@@ -177,4 +180,7 @@
       fCutCouple = theCoupleTable->GetMaterialCutsCouple( fMaterial, 
 					  curReg->GetProductionCuts() );
+      //G4cout << "Reg <" <<curReg->GetName() << ">  mat <" 
+      //	     << fMaterial->GetName() << ">  fCouple= " 
+      //	     << fCutCouple<<G4endl;
       if( fCutCouple ) {
 	fMaterialCutsCoupleVector.push_back(fCutCouple);
@@ -197,4 +203,9 @@
   }
 }
+
+//////////////////////////////////////////////////////////////////
+
+void G4PAIModel::InitialiseMe(const G4ParticleDefinition*) 
+{}
 
 //////////////////////////////////////////////////////////////////
@@ -393,5 +404,6 @@
   {
     //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*0.5 ;
     else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -435,5 +447,6 @@
   {
     //  if ( x1 == x2  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*0.5 ;
     else             dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -444,15 +457,20 @@
 //////////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIModel::ComputeDEDX(const G4MaterialCutsCouple* matCC,
-                                 const G4ParticleDefinition* p,
-                                       G4double kineticEnergy,
-                                       G4double cutEnergy)
+G4double G4PAIModel::ComputeDEDXPerVolume(const G4Material*,
+					  const G4ParticleDefinition* p,
+					  G4double kineticEnergy,
+					  G4double cutEnergy)
 {
   G4int iTkin,iPlace;
   size_t jMat;
+  
+  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
+  G4double cut = cutEnergy;
+
   G4double massRatio  = fMass/p->GetPDGMass();
   G4double scaledTkin = kineticEnergy*massRatio;
   G4double charge     = p->GetPDGCharge();
-  G4double charge2    = charge*charge, dEdx;
+  G4double charge2    = charge*charge;
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
@@ -470,6 +488,5 @@
   iPlace = iTkin - 1;
   if(iPlace < 0) iPlace = 0;
-  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ;  
-
+  G4double dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) );
   if( dEdx < 0.) dEdx = 0.;
   return dEdx;
@@ -478,17 +495,19 @@
 /////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIModel::CrossSection( const G4MaterialCutsCouple* matCC,
-                                   const G4ParticleDefinition* p,
-                                         G4double kineticEnergy,
-                                         G4double cutEnergy,
-                                         G4double maxEnergy  ) 
+G4double G4PAIModel::CrossSectionPerVolume( const G4Material*,
+					    const G4ParticleDefinition* p,
+					    G4double kineticEnergy,
+					    G4double cutEnergy,
+					    G4double maxEnergy  ) 
 {
   G4int iTkin,iPlace;
   size_t jMat;
-  G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  if(tmax <= cutEnergy) return 0.0;
   G4double massRatio  = fMass/p->GetPDGMass();
   G4double scaledTkin = kineticEnergy*massRatio;
   G4double charge     = p->GetPDGCharge();
   G4double charge2    = charge*charge, cross, cross1, cross2;
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
@@ -935,4 +954,27 @@
 }
 
+/////////////////////////////////////////////////////////////////////
+
+G4double G4PAIModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
+					 G4double kinEnergy) 
+{
+  G4double tmax = kinEnergy;
+  if(p == fElectron) tmax *= 0.5;
+  else if(p != fPositron) { 
+    G4double mass = p->GetPDGMass();
+    G4double ratio= electron_mass_c2/mass;
+    G4double gamma= kinEnergy/mass + 1.0;
+    tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
+                  (1. + 2.0*gamma*ratio + ratio*ratio);
+  }
+  return tmax;
+}
+
+///////////////////////////////////////////////////////////////
+
+void G4PAIModel::DefineForRegion(const G4Region* r) 
+{
+  fPAIRegionVector.push_back(r);
+}
 
 //
Index: trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4PAIPhotonModel.cc	(revision 1055)
@@ -24,4 +24,10 @@
 // ********************************************************************
 //
+// $Id: G4PAIPhotonModel.cc,v 1.22 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class 
 // File name:     G4PAIPhotonModel.cc
 //
@@ -160,8 +166,5 @@
   if(!fParticle) SetParticle(p);
 
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+  fParticleChange = GetParticleChangeForLoss();
 
   const G4ProductionCutsTable* theCoupleTable =
@@ -214,4 +217,9 @@
   }
 }
+
+//////////////////////////////////////////////////////////////////
+
+void G4PAIPhotonModel::InitialiseMe(const G4ParticleDefinition*)
+{}
 
 //////////////////////////////////////////////////////////////////
@@ -487,5 +495,6 @@
   {
     //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*0.5 ;
     else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -530,5 +539,6 @@
   {
     //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*0.5 ;
     else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -574,5 +584,6 @@
   {
     //  if ( x1 == x2  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dNdxCut = y1 + (y2 - y1)*0.5 ;
     else             dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -617,5 +628,6 @@
   {
     //  if ( x1 == x2  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
-    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    //    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ;
+    if ( std::abs(x1-x2) <= eV  ) dEdxCut = y1 + (y2 - y1)*0.5 ;
     else             dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ;      
   }
@@ -626,11 +638,15 @@
 //////////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIPhotonModel::ComputeDEDX(const G4MaterialCutsCouple* matCC,
-                                 const G4ParticleDefinition* p,
-                                       G4double kineticEnergy,
-                                       G4double cutEnergy)
+G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*,
+						const G4ParticleDefinition* p,
+						G4double kineticEnergy,
+						G4double cutEnergy)
 {
   G4int iTkin,iPlace;
   size_t jMat;
+
+  //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
+  G4double cut = cutEnergy;
+
   G4double particleMass = p->GetPDGMass();
   G4double scaledTkin   = kineticEnergy*proton_mass_c2/particleMass;
@@ -638,4 +654,5 @@
   G4double charge2      = charge*charge;
   G4double dEdx         = 0.;
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat )
@@ -653,5 +670,5 @@
   iPlace = iTkin - 1;
   if(iPlace < 0) iPlace = 0;
-  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cutEnergy) ) ;  
+  dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ;  
 
   if( dEdx < 0.) dEdx = 0.;
@@ -661,13 +678,14 @@
 /////////////////////////////////////////////////////////////////////////
 
-G4double G4PAIPhotonModel::CrossSection( const G4MaterialCutsCouple* matCC,
-                                   const G4ParticleDefinition* p,
-                                         G4double kineticEnergy,
-                                         G4double cutEnergy,
-                                         G4double maxEnergy  ) 
+G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*,
+						  const G4ParticleDefinition* p,
+						  G4double kineticEnergy,
+						  G4double cutEnergy,
+						  G4double maxEnergy  ) 
 {
   G4int iTkin,iPlace;
   size_t jMat, jMatCC;
-  G4double tmax = min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
+  if(cutEnergy >= tmax) return 0.0;
   G4double particleMass = p->GetPDGMass();
   G4double scaledTkin   = kineticEnergy*proton_mass_c2/particleMass;
@@ -675,4 +693,6 @@
   G4double charge2      = charge*charge, cross, cross1, cross2;
   G4double photon1, photon2, plasmon1, plasmon2;
+
+  const G4MaterialCutsCouple* matCC = CurrentCouple();
 
   const G4ProductionCutsTable* theCoupleTable=
@@ -1225,4 +1245,28 @@
 }
 
+/////////////////////////////////////////////////////////////////////
+
+G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
+                                                      G4double kinEnergy) 
+{
+  G4double tmax = kinEnergy;
+  if(p == fElectron) tmax *= 0.5;
+  else if(p != fPositron) { 
+    G4double mass = p->GetPDGMass();
+    G4double ratio= electron_mass_c2/mass;
+    G4double gamma= kinEnergy/mass + 1.0;
+    tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
+                  (1. + 2.0*gamma*ratio + ratio*ratio);
+  }
+  return tmax;
+}
+
+///////////////////////////////////////////////////////////////
+
+void G4PAIPhotonModel::DefineForRegion(const G4Region* r) 
+{
+  fPAIRegionVector.push_back(r);
+}
+
 
 //
Index: trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4PEEffectModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4PEEffectModel.cc,v 1.6 2007/05/22 17:34:36 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4PEEffectModel.cc,v 1.8 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -41,4 +41,6 @@
 //
 // 04.12.05 : SetProposedKineticEnergy(0.) for the killed photon (mma)
+// 20.02.09 : Added initialisation of deexcitation flag and method
+//            CrossSectionPerVolume instead of mfp (V.Ivanchenko)
 //
 // Class Description:
@@ -66,4 +68,5 @@
   theGamma    = G4Gamma::Gamma();
   theElectron = G4Electron::Electron();
+  fminimalEnergy = 1.0*eV;
 }
 
@@ -71,6 +74,5 @@
 
 G4PEEffectModel::~G4PEEffectModel()
-{
-}
+{}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -79,12 +81,45 @@
 				 const G4DataVector&)
 {
- if (isInitialized) return;
- if (pParticleChange)
-   fParticleChange =
-                  reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-  else
-   fParticleChange = new G4ParticleChangeForGamma();
-
- fminimalEnergy = 1.0*eV;
+  // always false before the run
+  SetDeexcitationFlag(false);
+
+  if (isInitialized) return;
+  fParticleChange = GetParticleChangeForGamma();
+  isInitialized = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4double G4PEEffectModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
+						     G4double energy,
+						     G4double Z, G4double,
+						     G4double, G4double)
+{
+  G4double* SandiaCof = G4SandiaTable::GetSandiaCofPerAtom((G4int)Z, energy);
+
+  G4double energy2 = energy*energy;
+  G4double energy3 = energy*energy2;
+  G4double energy4 = energy2*energy2;
+
+  return SandiaCof[0]/energy  + SandiaCof[1]/energy2 +
+    SandiaCof[2]/energy3 + SandiaCof[3]/energy4;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4PEEffectModel::CrossSectionPerVolume(const G4Material* material,
+						const G4ParticleDefinition*,
+						G4double energy,
+						G4double, G4double)
+{
+  G4double* SandiaCof = 
+    material->GetSandiaTable()->GetSandiaCofForMaterial(energy);
+				
+  G4double energy2 = energy*energy;
+  G4double energy3 = energy*energy2;
+  G4double energy4 = energy2*energy2;
+	  
+  return SandiaCof[0]/energy  + SandiaCof[1]/energy2 +
+    SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 
 }
 
Index: trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4PhotoElectricEffect.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4PhotoElectricEffect.cc,v 1.41 2008/10/16 14:12:32 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4PhotoElectricEffect.cc,v 1.42 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 //
@@ -89,4 +89,11 @@
 {}
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4PhotoElectricEffect::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Gamma::Gamma());
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
Index: trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4UniversalFluctuation.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UniversalFluctuation.cc,v 1.16 2008/10/22 16:04:33 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UniversalFluctuation.cc,v 1.22 2009/03/20 18:11:23 urban Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -57,5 +57,7 @@
 // 03-04-07 correction to get better width of eloss distr.(L.Urban)
 // 13-07-07 add protection for very small step or low-density material (VI)
-//          
+// 19-03-09 new width correction (does not depend on previous steps) (L.Urban)         
+// 20-03-09 modification in the width correction (L.Urban)
+//
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -84,7 +86,4 @@
 {
   lastMaterial = 0;
-  facwidth     = 1.000/keV;
-  oldloss = 0.;
-  samestep = 0.;
 }
 
@@ -120,8 +119,5 @@
   //
   if (meanLoss < minLoss)
-  { 
-    oldloss = meanLoss;
     return meanLoss;
-  }
 
   if(!particle) InitialiseMe(dp->GetDefinition());
@@ -179,4 +175,15 @@
     e0 = material->GetIonisation()->GetEnergy0fluct();
     lastMaterial = material;
+  
+    // modification of some model parameters
+    // (this part should go to materials later)
+    G4double p = 1.40;
+    f2Fluct *= p;
+    f1Fluct = 1.-f2Fluct;
+    G4double q = 1.00;
+    e2Fluct *= q;
+    e2LogFluct = log(e2Fluct);
+    e1LogFluct = (ipotLogFluct-f2Fluct*e2LogFluct)/f1Fluct;
+    e1Fluct = exp(e1LogFluct);
   }
 
@@ -185,15 +192,4 @@
 
   G4double a1 = 0. , a2 = 0., a3 = 0. ;
-
-  // correction to get better width even using stepmax
-  if(abs(meanLoss- oldloss) < 1.*eV)
-    samestep += 1;
-  else
-    samestep = 1.;
-  oldloss = meanLoss;
-  G4double width = 1.+samestep*facwidth*meanLoss;
-  if(width > 4.50) width = 4.50;
-  e1 = width*e1Fluct;
-  e2 = width*e2Fluct;
 
   // cut and material dependent rate 
@@ -206,6 +202,18 @@
       rate = 0.03+0.23*log(log(tmax/ipotFluct));
       G4double C = meanLoss*(1.-rate)/(w2-ipotLogFluct);
-      a1 = C*f1Fluct*(w2-e1LogFluct)/e1;
-      a2 = C*f2Fluct*(w2-e2LogFluct)/e2;
+      a1 = C*f1Fluct*(w2-e1LogFluct)/e1Fluct;
+      a2 = C*f2Fluct*(w2-e2LogFluct)/e2Fluct;
+      // correction in order to get better FWHM values
+      // ( scale parameters a1 and e1)
+      G4double width = 1.;
+      if(meanLoss > 10.*e1Fluct)
+      {
+        width = 3.1623/sqrt(meanLoss/e1Fluct);
+        if(width < a2/a1)
+        width = a2/a1;
+      } 
+      a1 *= width;
+      e1 = e1Fluct/width;
+      e2 = e2Fluct;
     }
   }
@@ -305,3 +313,16 @@
 }
 
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void 
+G4UniversalFluctuation::SetParticleAndCharge(const G4ParticleDefinition* part,
+					     G4double q2)
+{
+  if(part != particle) {
+    particle       = part;
+    particleMass   = part->GetPDGMass();
+  }
+  chargeSquare = q2;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel.cc	(revision 1055)
@@ -24,6 +24,7 @@
 // ********************************************************************
 //
-// $Id: G4UrbanMscModel.cc,v 1.86 2008/10/29 14:15:30 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+//
+// $Id: G4UrbanMscModel.cc,v 1.90 2009/04/29 13:30:22 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -36,118 +37,54 @@
 // Author:      Laszlo Urban
 //
-// Creation date: 03.03.2001
+// Creation date: 06.03.2008
 //
 // Modifications:
 //
-// 27-03-03 Move model part from G4MultipleScattering80 (V.Ivanchenko)
-// 23-05-03 important change in angle distribution for muons/hadrons
-//          the central part now is similar to the Highland parametrization +
-//          minor correction in angle sampling algorithm (for all particles)
-//          (L.Urban)
-// 30-05-03 misprint in SampleCosineTheta corrected(L.Urban)
-// 27-03-03 Rename (V.Ivanchenko)
-// 05-08-03 angle distribution has been modified (L.Urban)
-// 06-11-03 precision problems solved for high energy (PeV) particles
-//          change in the tail of the angular distribution
-//          highKinEnergy is set to 100 PeV (L.Urban) 
-//
-// 10-11-03 highKinEnergy is set back to 100 TeV, some tail tuning +
-//          cleaning (L.Urban) 
-// 26-11-03 correction in TrueStepLength : 
-//          trueLength <= currentRange (L.Urban) 
-// 01-03-04 signature changed in SampleCosineTheta,
-//          energy dependence calculations has been simplified,
-// 11-03-04 corrections in GeomPathLength,TrueStepLength,
-//          SampleCosineTheta
-// 23-04-04 true -> geom and geom -> true transformation has been
-//          rewritten, changes in the angular distribution (L.Urban)
-// 19-07-04 correction in SampleCosineTheta in order to avoid
-//          num. precision problems at high energy/small step(L.Urban) 
-// 17-08-04 changes in the angle distribution (slightly modified
-//          Highland formula for the width of the central part,
-//          changes in the numerical values of some other parameters)
-//          ---> approximately step independent distribution (L.Urban)
-// 21-09-04 change in the tail of the angular distribution (L.Urban)
-//
-// 03-11-04 precision problem for very high energy ions and small stepsize
-//          solved in SampleCosineTheta (L.Urban).
-// 15-04-05 optimize internal interface 
-//          add SampleSecondaries method (V.Ivanchenko)
-// 11-08-05 computation of lateral correlation added (L.Urban)
-// 02-10-05 nuclear size correction computation removed, the correction
-//          included in the (theoretical) tabulated values (L.Urban)
-// 17-01-06 computation of tail changed in SampleCosineTheta (l.Urban)
-// 16-02-06 code cleaning + revised 'z' sampling (L.Urban)
-// 17-02-06 Save table of transport cross sections not mfp (V.Ivanchenko)
-// 07-03-06 Create G4UrbanMscModel and move there step limit
-//           calculation (V.Ivanchenko)
-// 23-03-06 Bugfix in SampleCosineTheta method (L.Urban)
-// 10-05-06 SetMscStepLimitation at initialisation (V.Ivantchenko)
-// 11-05-06 name of data member safety changed to presafety, some new data
-//          members added (frscaling1,frscaling2,tlimitminfix,nstepmax)
-//          changes in ComputeTruePathLengthLimit,SampleCosineTheta (L.Urban)
-// 17-05-06 parameters of theta0 in SampleCosineTheta changed
-//          c_highland  13.6*MeV ---> 13.26*MeV,
-//          corr_highland  0.555 ---> 0.54,
-//          value of data member geommin changed (5 nm -> 1 nm) (L.Urban)
-// 13-10-06 data member factail removed, data member tkinlimit changed
-//          to lambdalimit,
-//          new data members tgeom,tnow,skin,skindepth,Zeff,geomlimit
-//          G4double GeomLimit(const G4Track& track) changed to
-//              void GeomLimit(const G4Track& track) 
-//        - important changes in ComputeTruePathLengthLimit: 
-//          possibility to have very small step(s) with single scattering 
-//          before boundary crossing (with skin > 0)
-//        - changes in SampleCosineTheta :
-//          single scattering if step <= stepmin, parameter theta0 
-//          slightly modified, tail modified (L.Urban)
-// 20-10-06 parameter theta0 now computed in the (public)
-//          function ComputeTheta0,
-//          single scattering modified allowing not small
-//          angles as well (L.Urban)
-// 23-10-06 correction in SampleSecondaries, now safety update
-//          computed in a simpler/faster way (L.Urban)
-// 06-11-06 corrections in ComputeTruePathLengthLimit, results are
-//          more stable in calorimeters (L.Urban)
-// 07-11-06 fix in GeomPathLength and SampleCosineTheta (L.Urban)
-// 15-11-06 bugfix in SampleCosineTheta (L.Urban)
-// 20-11-06 bugfix in single scattering part of SampleCosineTheta,
-//          single scattering just before boundary crossing now (L.Urban)
-// 04-12-06 fix in ComputeTruePathLengthLimit (L.Urban)
-// 17-01-07 remove LocatePoint from GeomLimit method (V.Ivanchenko)
-// 19-01-07 fix of true < geom problem (L.Urban)
-// 25-01-07 add protections from NaN vaues and for zero geometry step (VI)
-// 31-01-07 correction in SampleCosineTheta: screening parameter 
-//          corrected in single/plural scattering +
-//          code cleaning (L.Urban)
-// 01-02-07 restore logic inside ComputeTrueStepLength (V.Ivanchenko)
-// 06-02-07 Move SetMscStepLimitation method into the source, add there 
-//          reinitialisation of some private members, add protection inside 
-//          SampleDisplacement(VI)
-// 07-02-07 fix single scattering for heavy particles, now skin=1 can be used
-//          for heavy particles as well (L.Urban)
-// 08-02-07 randomization of tlimit removed (L.Urban)
-// 11-02-07 modified stepping algorithm for skin=0
-// 15-02-07 new data member: smallstep, small steps with single scattering
-//          before + after boundary for skin > 1
-// 23-02-07 use tPathLength inside ComputeStep instead of geombig
-// 24-02-07 step reduction before boundary for 'small' geomlimit only 
-// 03-03-07 single scattering around boundaries only (L.Urban)
-// 07-03-07 bugfix in ComputeTruePathLengthLimit (for skin > 0.) (L.Urban)
-// 10-04-07 optimize logic of ComputeTruePathLengthLimit, remove
-//          unused members, use unique G4SafetyHelper (V.Ivanchenko)
-// 01-05-07 optimization for skin > 0 (L.Urban)
-// 05-07-07 modified model functions in SampleCosineTheta (L.Urban)
-// 06-07-07 theta0 is not the same for e-/e+ as for heavy particles (L.Urban)
-// 02-08-07 compare safety not with 0. but with tlimitminfix (V.Ivanchenko)
-// 09-08-07 tail of angular distribution has been modified (L.Urban)
-// 22-10-07 - corr. in ComputeGeomPathLength in order to get better low
-//          energy behaviour for heavy particles,
-//          - theta0 is slightly modified,
-//          - some old inconsistency/bug is cured in SampleCosineTheta,
-//          now 0 <= prob <= 1 in any case (L.Urban) 
-// 26-10-07 different correction parameters for e/mu/hadrons in ComputeTheta0
-// 30-11-07 fix in ComputeTheta0 (L.Urban)
-//
+// 06-03-2008 starting point : G4UrbanMscModel2 = G4UrbanMscModel 9.1 ref 02
+//
+// 13-03-08  Bug in SampleScattering (which caused lateral asymmetry) fixed
+//           (L.Urban)
+//
+// 14-03-08  Simplification of step limitation in ComputeTruePathLengthLimit,
+//           + tlimitmin is the same for UseDistancetoBoundary and
+//           UseSafety (L.Urban)           
+//
+// 16-03-08  Reorganization of SampleCosineTheta + new method SimpleScattering
+//           SimpleScattering is used if the relative energy loss is too big
+//           or theta0 is too big (see data members rellossmax, theta0max)
+//           (L.Urban)          
+//
+// 17-03-08  tuning of the correction factor in ComputeTheta0 (L.Urban)
+//
+// 19-03-08  exponent c of the 'tail' model function is not equal to 2 any more,
+//           value of c has been extracted from some e- scattering data (L.Urban) 
+//
+// 24-03-08  Step limitation in ComputeTruePathLengthLimit has been
+//           simplified further + some data members have been removed (L.Urban)
+//
+// 24-07-08  central part of scattering angle (theta0) has been tuned
+//           tail of the scattering angle distribution has been tuned
+//           using some e- and proton scattering data
+//
+// 05-08-08  bugfix in ComputeTruePathLengthLimit (L.Urban)
+//
+// 09-10-08  theta0 and tail have been retuned using some e-,mu,proton
+//           scattering data (L.Urban)
+//           + single scattering without path length correction for
+//           small steps (t < tlimitmin, for UseDistanceToBoundary only)
+//
+// 15-10-08  Moliere-Bethe screening in the single scattering part(L.Urban)          
+//
+// 17-10-08  stepping similar to that in model (9.1) for UseSafety case
+//           for e+/e- in order to speed up the code for calorimeters
+//
+// 23-10-08  bugfix in the screeningparameter of the single scattering part,
+//           some technical change in order to speed up the code (UpdateCache)
+//
+// 27-10-08  bugfix in ComputeTruePathLengthLimit (affects UseDistanceToBoundary
+//           stepping type only) (L.Urban)          
+//
+// 28-04-09  move G4UrbanMscModel2 from the g49.2 to G4UrbanMscModel. 
+//           now it is frozen (V.Ivanchenk0)
 
 // Class Description:
@@ -157,4 +94,6 @@
 
 // -------------------------------------------------------------------
+// In its present form the model can be  used for simulation 
+//   of the e-/e+, muon and charged hadron multiple scattering
 //
 
@@ -166,11 +105,9 @@
 #include "Randomize.hh"
 #include "G4Electron.hh"
-
 #include "G4LossTableManager.hh"
 #include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
 
 #include "G4Poisson.hh"
+#include "globals.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -182,7 +119,8 @@
     isInitialized(false)
 {
-  masslimite  = 0.6*MeV;
-  masslimitmu = 110.*MeV;
-
+  masslimite    = 0.6*MeV;
+  lambdalimit   = 1.*mm;
+  fr            = 0.02;
+  facsafety     = 0.3;
   taubig        = 8.0;
   tausmall      = 1.e-16;
@@ -193,14 +131,33 @@
   smallstep     = 1.e10;
   currentRange  = 0. ;
-  frscaling2    = 0.25;
-  frscaling1    = 1.-frscaling2;
+  rangeinit     = 0.;
   tlimit        = 1.e10*mm;
   tlimitmin     = 10.*tlimitminfix;            
-  nstepmax      = 25.;
+  tgeom         = 1.e50*mm;
   geombig       = 1.e50*mm;
   geommin       = 1.e-3*mm;
   geomlimit     = geombig;
   presafety     = 0.*mm;
+                          
+  y             = 0.;
+
+  Zold          = 0.;
   Zeff          = 1.;
+  Z2            = 1.;                
+  Z23           = 1.;                    
+  lnZ           = 0.;
+  coeffth1      = 0.;
+  coeffth2      = 0.;
+  coeffc1       = 0.;
+  coeffc2       = 0.;
+  scr1ini       = fine_structure_const*fine_structure_const*
+                  electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
+  scr2ini       = 3.76*fine_structure_const*fine_structure_const;
+  scr1          = 0.;
+  scr2          = 0.;
+
+  theta0max     = pi/6.;
+  rellossmax    = 0.50;
+  third         = 1./3.;
   particle      = 0;
   theManager    = G4LossTableManager::Instance(); 
@@ -218,20 +175,13 @@
 
 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p,
-				 const G4DataVector&)
+				  const G4DataVector&)
 {
   skindepth = skin*stepmin;
   if(isInitialized) return;
-
   // set values of some data members
   SetParticle(p);
 
-  if (pParticleChange)
-   fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-  else
-   fParticleChange = new G4ParticleChangeForMSC();
-
-  safetyHelper = G4TransportationManager::GetTransportationManager()
-    ->GetSafetyHelper();
-  safetyHelper->InitialiseHelper();
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
 
   isInitialized = true;
@@ -385,12 +335,10 @@
   if(mass > electron_mass_c2)
   {
-    G4double TAU = KineticEnergy/mass ;
-    G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
-    G4double w = c-2. ;
-    G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
-    eKineticEnergy = electron_mass_c2*tau ;
-  }
-
-  G4double ChargeSquare = charge*charge;
+     G4double TAU = KineticEnergy/mass ;
+     G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
+     G4double w = c-2. ;
+     G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
+     eKineticEnergy = electron_mass_c2*tau ;
+  }
 
   G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
@@ -490,7 +438,8 @@
 {
   tPathLength = currentMinimalStep;
-  const G4DynamicParticle* dp = track.GetDynamicParticle();
   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
   G4StepStatus stepStatus = sp->GetStepStatus();
+
+  const G4DynamicParticle* dp = track.GetDynamicParticle();
 
   if(stepStatus == fUndefined) {
@@ -516,7 +465,7 @@
   presafety = sp->GetSafety();
 
-  //  G4cout << "G4UrbanMscModel::ComputeTruePathLengthLimit tPathLength= " 
-  //	 <<tPathLength<<" safety= " << presafety
-  //     << " range= " <<currentRange<<G4endl;
+  //    G4cout << "G4UrbanMscModel::ComputeTruePathLengthLimit tPathLength= " 
+  //	   <<tPathLength<<" safety= " << presafety
+  //	   << " range= " <<currentRange<<G4endl;
 
   // far from geometry boundary
@@ -532,8 +481,8 @@
     {
       //compute geomlimit and presafety 
-      GeomLimit(track);
-
-      // is far from boundary
-      if(currentRange <= presafety)
+      G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
+
+      // is it far from boundary ?
+      if(currentRange < presafety)
 	{
 	  inside = true;
@@ -546,37 +495,9 @@
       if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
 	{
+          rangeinit = currentRange;
 	  if(stepStatus == fUndefined) smallstep = 1.e10;
 	  else  smallstep = 1.;
 
-	  // facrange scaling in lambda 
-	  // not so strong step restriction above lambdalimit
-	  G4double facr = facrange;
-	  if(lambda0 > lambdalimit)
-	    facr *= frscaling1+frscaling2*lambda0/lambdalimit;
-
-	  // constraint from the physics
-	  if (currentRange > lambda0) tlimit = facr*currentRange;
-	  else                        tlimit = facr*lambda0;
-
-          if(tlimit > currentRange) tlimit = currentRange;
-
-	  //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*lambda0;
-	  skindepth = skin*stepmin;
-
-	  //define tlimitmin
-          tlimitmin = 10.*stepmin;
-	  if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
-
-	  //lower limit for tlimit
-	  if(tlimit < tlimitmin) tlimit = tlimitmin;
-
-	  // constraint from the geometry (if tlimit above is too big)
-	  G4double tgeom = geombig; 
-
+	  // constraint from the geometry 
 	  if((geomlimit < geombig) && (geomlimit > geommin))
 	    {
@@ -585,11 +506,31 @@
 	      else
 	        tgeom = 2.*geomlimit/facgeom;
-
-	      if(tlimit > tgeom) tlimit = tgeom;
 	    }
-	}
-
-      //if track starts far from boundaries increase tlimit!
-      if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ;
+            else
+              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*lambda0;
+          skindepth = skin*stepmin;
+
+          //define tlimitmin
+          tlimitmin = 10.*stepmin;
+          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+
+        }
+
+      //step limit 
+      tlimit = facrange*rangeinit;              
+      if(tlimit < facsafety*presafety)
+        tlimit = facsafety*presafety; 
+
+      //lower limit for tlimit
+      if(tlimit < tlimitmin) tlimit = tlimitmin;
+
+      if(tlimit > tgeom) tlimit = tgeom;
 
       //  G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
@@ -597,15 +538,12 @@
 
       // shortcut
-      if((tPathLength < tlimit) && (tPathLength < presafety))
+      if((tPathLength < tlimit) && (tPathLength < presafety) &&
+         (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
 	return tPathLength;   
-
-      G4double tnow = tlimit;
-      // optimization ...
-      if(geomlimit < geombig) tnow = max(tlimit,facsafety*geomlimit);
 
       // step reduction near to boundary
       if(smallstep < skin)
 	{
-	  tnow = stepmin;
+	  tlimit = stepmin;
 	  insideskin = true;
 	}
@@ -614,17 +552,18 @@
 	  if(geomlimit > skindepth)
 	    {
-	      if(tnow > geomlimit-0.999*skindepth)
-		tnow = geomlimit-0.999*skindepth;
+	      if(tlimit > geomlimit-0.999*skindepth)
+		tlimit = geomlimit-0.999*skindepth;
 	    }
 	  else
 	    {
 	      insideskin = true;
-	      if(tnow > stepmin) tnow = stepmin;
+	      if(tlimit > stepmin) tlimit = stepmin;
 	    }
 	}
 
-      if(tnow < stepmin) tnow = stepmin;
-
-      if(tPathLength > tnow) tPathLength = tnow ; 
+      if(tlimit < stepmin) tlimit = stepmin;
+
+      if(tPathLength > tlimit) tPathLength = tlimit  ; 
+
     }
     // for 'normal' simulation with or without magnetic field 
@@ -635,5 +574,5 @@
       // i.e. when it is needed for optimization purposes
       if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
-        presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
+	presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
 
       // is far from boundary
@@ -645,22 +584,30 @@
 
       if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
-	{ 
-	  // facrange scaling in lambda 
-	  // not so strong step restriction above lambdalimit
-	  G4double facr = facrange;
-	  if(lambda0 > lambdalimit)
-	    facr *= frscaling1+frscaling2*lambda0/lambdalimit;
-
-	  // constraint from the physics
-	  if (currentRange > lambda0) tlimit = facr*currentRange;
-	  else                        tlimit = facr*lambda0;
-
-	  //lower limit for tlimit
-	  tlimitmin = std::max(tlimitminfix,lambda0/nstepmax);
-	  if(tlimit < tlimitmin) tlimit = tlimitmin;
-	}
-
-      //if track starts far from boundaries increase tlimit!
-      if(tlimit < facsafety*presafety) tlimit = facsafety*presafety ;
+      {
+        rangeinit = currentRange;
+        fr = facrange;
+        // 9.1 like stepping for e+/e- only (not for muons,hadrons)
+        if(mass < masslimite) 
+        {
+          if(lambda0 > currentRange)
+            rangeinit = lambda0;
+          if(lambda0 > lambdalimit)
+            fr *= 0.75+0.25*lambda0/lambdalimit;
+        }
+
+        //lower limit for tlimit
+        G4double rat = currentKinEnergy/MeV ;
+        rat = 1.e-3/(rat*(10.+rat)) ;
+        tlimitmin = 10.*lambda0*rat;
+        if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+      }
+      //step limit
+      tlimit = fr*rangeinit;               
+
+      if(tlimit < facsafety*presafety)
+        tlimit = facsafety*presafety;
+
+      //lower limit for tlimit
+      if(tlimit < tlimitmin) tlimit = tlimitmin;
 
       if(tPathLength > tlimit) tPathLength = tlimit;
@@ -679,30 +626,7 @@
 	}
     }
-  //  G4cout << "tPathLength= " << tPathLength << "  geomlimit= " << geomlimit 
+  // G4cout << "tPathLength= " << tPathLength << "  geomlimit= " << geomlimit 
   //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
-
   return tPathLength ;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4UrbanMscModel::GeomLimit(const G4Track&  track)
-{
-  geomlimit = geombig;
-
-  // no geomlimit for the World volume
-  if((track.GetVolume() != 0) &&
-     (track.GetVolume() != safetyHelper->GetWorldVolume()))  
-  {
-    G4double cstep = currentRange;
-
-    geomlimit = safetyHelper->CheckNextStep(
-                  track.GetStep()->GetPreStepPoint()->GetPosition(),
-                  track.GetMomentumDirection(),
-                  cstep,
-                  presafety);
-    //    G4cout << "!!!G4UrbanMscModel::GeomLimit presafety= " << presafety
-    //	   << " limit= " << geomlimit << G4endl;
-  }  
 }
 
@@ -761,5 +685,5 @@
   if(samplez)
   {
-    const G4double  ztmax = 0.99, onethird = 1./3. ;
+    const G4double  ztmax = 0.99 ;
     G4double zt = zmean/tPathLength ;
 
@@ -767,5 +691,5 @@
     {
       G4double u,cz1;
-      if(zt >= onethird)
+      if(zt >= third)
       {
         G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
@@ -835,14 +759,12 @@
                          KineticEnergy*(KineticEnergy+2.*mass)/
                       ((currentKinEnergy+mass)*(KineticEnergy+mass)));
-  G4double y = trueStepLength/currentRadLength;
+  y = trueStepLength/currentRadLength;
   G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
   y = log(y);
-  if(mass < masslimite)                 
-    theta0 *= (1.+0.051*y);
-  else if(mass < masslimitmu)
-    theta0 *= (1.+0.044*y);
-  else
-    theta0 *= (1.+0.038*y);
-    
+  // correction factor from e-/proton scattering data
+  G4double corr = coeffth1+coeffth2*y;                
+  if(y < -6.5) corr -= 0.011*(6.5+y);
+  theta0 *= corr ;                                               
+
   return theta0;
 }
@@ -854,10 +776,12 @@
 {
   G4double kineticEnergy = dynParticle->GetKineticEnergy();
+
   if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) ||
-     (tPathLength/tausmall < lambda0) ) return;
+     (tPathLength/tausmall < lambda0)) return;
 
   G4double cth  = SampleCosineTheta(tPathLength,kineticEnergy);
+
   // protection against 'bad' cth values
-  if(abs(cth) > 1.) return;
+  if(std::abs(cth) > 1.) return;
 
   G4double sth  = sqrt((1.0 - cth)*(1.0 + cth));
@@ -874,5 +798,5 @@
 
     G4double r = SampleDisplacement();
-/*
+    /*
     G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
 	   << " sinTheta= " << sth << " r(mm)= " << r
@@ -880,5 +804,5 @@
            << " geomStep(mm)= " << zPathLength
            << G4endl;
-*/
+    */
     if(r > 0.)
       {
@@ -889,11 +813,14 @@
         // compute it from the lateral correlation
         G4double Phi = 0.;
-        if(std::abs(r*sth) < latcorr) {
+        if(std::abs(r*sth) < latcorr)
           Phi  = twopi*G4UniformRand();
-        } else {
+        else
+        {
           G4double psi = std::acos(latcorr/(r*sth));
-          if(G4UniformRand() < 0.5) Phi = phi+psi;
-          else                      Phi = phi-psi;
-	}
+          if(G4UniformRand() < 0.5)
+            Phi = phi+psi;
+          else
+            Phi = phi-psi;
+        }
 
         dirx = std::cos(Phi);
@@ -903,39 +830,6 @@
         latDirection.rotateUz(oldDirection);
 
-        G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-        G4double fac = 1.;
-        if(r >  safety) {
-          //  ******* so safety is computed at boundary too ************
-	  G4double newsafety = safetyHelper->ComputeSafety(Position);
-          if(r > newsafety)
-            fac = newsafety/r ;
-        }  
-
-        if(fac > 0.)
-        {
-          // compute new endpoint of the Step
-          G4ThreeVector newPosition = Position+fac*r*latDirection;
-
-	  // definitely not on boundary
-	  if(1. == fac) {
-	    safetyHelper->ReLocateWithinVolume(newPosition);
-	    
-	  } else {
-            // check safety after displacement
-	    G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	    // displacement to boundary
-	    // if(postsafety < tlimitminfix) {
-            if(postsafety <= 0.0) {
-	      safetyHelper->Locate(newPosition, newDirection);
-
-	    // not on the boundary
-            } else { 
-	      safetyHelper->ReLocateWithinVolume(newPosition);
-	    }
-	  }
-          fParticleChange->ProposePosition(newPosition);
-        } 
-     }
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }
   }
 }
@@ -951,4 +845,7 @@
   Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
          couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
+
+  if(Zold != Zeff)  
+    UpdateCache();
 
   if(insideskin)
@@ -960,10 +857,13 @@
     if(n > 0)
     {
-      G4double tm = KineticEnergy/electron_mass_c2;
-      // ascr - screening parameter
-      G4double ascr = exp(log(Zeff)/3.)/(137.*sqrt(tm*(tm+2.)));
-      G4double ascr1 = 1.+0.5*ascr*ascr;
+      //screening (Moliere-Bethe)
+      G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
+      G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
+      G4double ascr = scr1/mom2;
+      ascr *= 1.13+scr2/beta2;
+      G4double ascr1 = 1.+2.*ascr;
       G4double bp1=ascr1+1.;
       G4double bm1=ascr1-1.;
+
       // single scattering from screened Rutherford x-section
       G4double ct,st,phi;
@@ -1001,75 +901,92 @@
     else if (tau >= tausmall)
     {
-      G4double b=2.,bp1=3.,bm1=1.;
-      G4double prob = 0. ;
+      G4double xsi = 3.0;  
+      G4double x0 = 1.;
       G4double a = 1., ea = 0., eaa = 1.;
+      G4double b=2.,b1=3.,bx=1.,eb1=3.,ebx=1.;
+      G4double prob = 1. , qprob = 1. ;
       G4double xmean1 = 1., xmean2 = 0.;
-                                                      
+      G4double xmeanth = exp(-tau);
+      G4double x2meanth = (1.+2.*exp(-2.5*tau))/3.;
+
+      G4double relloss = 1.-KineticEnergy/currentKinEnergy;
+      if(relloss > rellossmax) 
+        return SimpleScattering(xmeanth,x2meanth);
+
       G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
 
-      // protexction for very small angles
+      // protection for very small angles
       if(theta0 < tausmall) return cth;
-
+    
+      if(theta0 > theta0max)
+        return SimpleScattering(xmeanth,x2meanth);
       G4double sth = sin(0.5*theta0);
       a = 0.25/(sth*sth);
 
-      ea = exp(-2.*a);
-      eaa = 1.-ea ; 
-      xmean1 = (1.+ea)/eaa-1./a;
-
-      G4double xmeanth = exp(-tau);
-      b = 1./xmeanth ;                               
-      bp1 = b+1.;
-      bm1 = b-1.;
-      // protection 
-      if(bm1 > 0.)
-        xmean2 = b-0.5*bp1*bm1*log(bp1/bm1);
-      else
+      ea = exp(-xsi);
+      eaa = 1.-ea ;
+      xmean1 = 1.-(1.-(1.+xsi)*ea)/(a*eaa);
+      x0 = 1.-xsi/a;
+
+      if(xmean1 <= 0.999*xmeanth)
+        return SimpleScattering(xmeanth,x2meanth);
+
+      // from MUSCAT H,Be,Fe data 
+      G4double c = coeffc1;                         
+      if(y > -13.5) 
+        c += coeffc2*exp(3.*log(y+13.5));
+
+      if(abs(c-3.) < 0.001)  c = 3.001;      
+      if(abs(c-2.) < 0.001)  c = 2.001;      
+
+      G4double c1 = c-1.;
+
+      //from continuity of derivatives
+      b = 1.+(c-xsi)/a;
+
+      b1 = b+1.;
+      bx = c/a;
+      eb1 = exp(c1*log(b1));
+      ebx = exp(c1*log(bx));
+
+      xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx);
+      
+      G4double f1x0 = a*ea/eaa;
+      G4double f2x0 = c1*eb1/(bx*(eb1-ebx));
+      prob = f2x0/(f1x0+f2x0);
+
+      qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
+
+      // sampling of costheta
+      if(G4UniformRand() < qprob)
       {
-        b = 1.+tau;
-        bp1 = 2.+tau;
-        bm1 = tau;
-        xmean2 = 1.+tau*(1.-log(2./tau));
-      }
-
-      if((xmean1 >= xmeanth) && (xmean2 <= xmeanth))
-      {
-        //normal case
-        prob = (xmeanth-xmean2)/(xmean1-xmean2);
+        if(G4UniformRand() < prob)
+          cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
+        else
+          cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
       }
       else
-      {
-        // x1 < xth ( x2 < xth automatically if b = 1/xth)
-        // correct a (xmean1)
-        if((xmeanth-xmean1)/xmeanth < 1.e-5)
-        {   
-          // xmean1 is small probably due to precision problems
-          xmean1 = 0.50*(1.+xmeanth);
-          prob = (xmeanth-xmean2)/(xmean1-xmean2);
-        }
-        else
-        {
-          //  correct a in order to have x1=xth
-          G4int i=0, imax=10;
-          do
-          {
-            a = 1./(1.-xmeanth+2.*ea/eaa);
-            ea = exp(-2.*a);
-            eaa = 1.-ea;
-            xmean1 = (1.+ea)/eaa-1./a;
-            i += 1;
-          } while ((std::abs((xmeanth-xmean1)/xmeanth) > 0.05) && (i < imax));
-	  prob = 1.;          
-        }
-      }
-
-      // sampling of costheta
-      if (G4UniformRand() < prob)
-         cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
-      else
-         cth = b-bp1*bm1/(bm1+2.*G4UniformRand()) ;
+        cth = -1.+2.*G4UniformRand();
     }
   }  
   return cth ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel::SimpleScattering(G4double xmeanth,G4double x2meanth)
+{
+  // 'large angle scattering'
+  // 2 model functions with correct xmean and x2mean
+  G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
+  G4double prob = (a+2.)*xmeanth/a;
+
+  // sampling
+  G4double cth = 1.;
+  if(G4UniformRand() < prob)
+    cth = -1.+2.*exp(log(G4UniformRand())/(a+1.));
+  else
+    cth = -1.+2.*G4UniformRand();
+  return cth;
 }
 
@@ -1137,11 +1054,2 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4UrbanMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
-					const G4MaterialCutsCouple*,
-					const G4DynamicParticle*,
-					G4double,
-					G4double)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel2.cc	(revision 1055)
@@ -25,6 +25,6 @@
 //
 //
-// $Id: G4UrbanMscModel2.cc,v 1.18 2008/12/18 13:01:36 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UrbanMscModel2.cc,v 1.26 2009/05/19 06:26:10 urban Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -104,6 +104,4 @@
 #include "G4LossTableManager.hh"
 #include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
 
 #include "G4Poisson.hh"
@@ -181,12 +179,6 @@
   SetParticle(p);
 
-  if (pParticleChange)
-   fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-  else
-   fParticleChange = new G4ParticleChangeForMSC();
-
-  safetyHelper = G4TransportationManager::GetTransportationManager()
-    ->GetSafetyHelper();
-  safetyHelper->InitialiseHelper();
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
 
   isInitialized = true;
@@ -486,5 +478,5 @@
     {
       //compute geomlimit and presafety 
-      GeomLimit(track);
+      G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
 
       // is it far from boundary ?
@@ -499,19 +491,8 @@
 
       if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
-	{
+        {
           rangeinit = currentRange;
-	  if(stepStatus == fUndefined) smallstep = 1.e10;
-	  else  smallstep = 1.;
-
-	  // constraint from the geometry 
-	  if((geomlimit < geombig) && (geomlimit > geommin))
-	    {
-	      if(stepStatus == fGeomBoundary)  
-	        tgeom = geomlimit/facgeom;
-	      else
-	        tgeom = 2.*geomlimit/facgeom;
-	    }
-            else
-              tgeom = geombig;
+          if(stepStatus == fUndefined) smallstep = 1.e10;
+          else  smallstep = 1.;
 
           //define stepmin here (it depends on lambda!)
@@ -522,10 +503,25 @@
           stepmin = rat*lambda0;
           skindepth = skin*stepmin;
-
           //define tlimitmin
           tlimitmin = 10.*stepmin;
           if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
 
+          // constraint from the geometry
+          if((geomlimit < geombig) && (geomlimit > geommin))
+            {
+              // geomlimit is a geometrical step length
+              // transform it to true path length (estimation)
+              if((1.-geomlimit/lambda0) > 0.)
+                geomlimit = -lambda0*log(1.-geomlimit/lambda0)+tlimitmin ;
+
+              if(stepStatus == fGeomBoundary)
+                tgeom = geomlimit/facgeom;
+              else
+                tgeom = 2.*geomlimit/facgeom;
+            }
+            else
+              tgeom = geombig;
         }
+
 
       //step limit 
@@ -579,5 +575,5 @@
       // i.e. when it is needed for optimization purposes
       if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
-        presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
+	presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
 
       // is far from boundary
@@ -634,29 +630,4 @@
   //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
   return tPathLength ;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4UrbanMscModel2::GeomLimit(const G4Track&  track)
-{
-  geomlimit = geombig;
-
-  // no geomlimit for the World volume
-  if((track.GetVolume() != 0) &&
-     (track.GetVolume() != safetyHelper->GetWorldVolume()))  
-  {
-    G4double cstep = currentRange;
-
-    geomlimit = safetyHelper->CheckNextStep(
-                  track.GetStep()->GetPreStepPoint()->GetPosition(),
-                  track.GetMomentumDirection(),
-                  cstep,
-                  presafety);
-    /*
-    G4cout << "!!!G4UrbanMscModel2::GeomLimit presafety= " << presafety
-	   << " r= " << currentRange
-    	   << " limit= " << geomlimit << G4endl;
-    */
-  }  
 }
 
@@ -792,7 +763,7 @@
   G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
   y = log(y);
-  // correction factor from e-/proton scattering data
+  // correction factor from e- scattering data
   G4double corr = coeffth1+coeffth2*y;                
-  if(y < -6.5) corr -= 0.011*(6.5+y);
+
   theta0 *= corr ;                                               
 
@@ -860,38 +831,6 @@
         latDirection.rotateUz(oldDirection);
 
-        G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-        G4double fac = 1.;
-        if(r >  safety) {
-          //  ******* so safety is computed at boundary too ************
-	  G4double newsafety = safetyHelper->ComputeSafety(Position);
-          if(r > newsafety)
-            fac = newsafety/r ;
-        }  
-
-        if(fac > 0.)
-        {
-          // compute new endpoint of the Step
-          G4ThreeVector newPosition = Position+fac*r*latDirection;
-
-	  // definitely not on boundary
-	  if(1. == fac) {
-	    safetyHelper->ReLocateWithinVolume(newPosition);
-	    
-	  } else {
-            // check safety after displacement
-	    G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	    // displacement to boundary
-            if(postsafety <= 0.0) {
-	      safetyHelper->Locate(newPosition, newDirection);
-
-	    // not on the boundary
-            } else {
-	      safetyHelper->ReLocateWithinVolume(newPosition);
-	    }
-	  }
-          fParticleChange->ProposePosition(newPosition);
-        } 
-     }
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }
   }
 }
@@ -963,5 +902,5 @@
     else if (tau >= tausmall)
     {
-      G4double xsi = 3.0;  
+      G4double xsi = 3.;
       G4double x0 = 1.;
       G4double a = 1., ea = 0., eaa = 1.;
@@ -994,11 +933,10 @@
         return SimpleScattering(xmeanth,x2meanth);
 
-      // from MUSCAT H,Be,Fe data 
-      G4double c = coeffc1;                         
-      if(y > -13.5) 
-        c += coeffc2*exp(3.*log(y+13.5));
+      // from e- and muon scattering data                    
+      G4double c = coeffc1+coeffc2*y; ;                         
 
       if(abs(c-3.) < 0.001)  c = 3.001;      
       if(abs(c-2.) < 0.001)  c = 2.001;      
+      if(abs(c-1.) < 0.001)  c = 1.001;      
 
       G4double c1 = c-1.;
@@ -1116,11 +1054,2 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4UrbanMscModel2::SampleSecondaries(std::vector<G4DynamicParticle*>*,
-					 const G4MaterialCutsCouple*,
-					 const G4DynamicParticle*,
-					 G4double,
-					 G4double)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel90.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4UrbanMscModel90.cc,v 1.10 2008/10/29 14:15:30 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4UrbanMscModel90.cc,v 1.13 2009/04/10 18:10:58 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -60,6 +60,4 @@
 #include "G4LossTableManager.hh"
 #include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
 
 #include "G4Poisson.hh"
@@ -113,12 +111,6 @@
   SetParticle(p);
 
-  if (pParticleChange) {
-   fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-  } else {
-   fParticleChange = new G4ParticleChangeForMSC();
-  }
-  safetyHelper = G4TransportationManager::GetTransportationManager()
-    ->GetSafetyHelper();
-  safetyHelper->InitialiseHelper();
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
 
   isInitialized = true;
@@ -419,5 +411,5 @@
     {
       //compute geomlimit and presafety 
-      GeomLimit(track);
+      G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
    
       // is far from boundary
@@ -522,5 +514,5 @@
       // i.e. when it is needed for optimization purposes
       if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
-        presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
+	presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
 
       // is far from boundary
@@ -569,25 +561,4 @@
   //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
   return tPathLength ;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4UrbanMscModel90::GeomLimit(const G4Track&  track)
-{
-  geomlimit = geombig;
-
-  // no geomlimit for the World volume
-  if((track.GetVolume() != 0) &&
-     (track.GetVolume() != safetyHelper->GetWorldVolume()))  
-  {
-    G4double cstep = currentRange;
-    geomlimit = safetyHelper->CheckNextStep(
-                  track.GetStep()->GetPreStepPoint()->GetPosition(),
-                  track.GetMomentumDirection(),
-                  cstep,
-                  presafety);
-    //    G4cout << "!!!G4UrbanMscModel90::GeomLimit presafety= " << presafety
-    //	   << " limit= " << geomlimit << G4endl;
-  }  
 }
 
@@ -788,42 +759,6 @@
         G4ThreeVector latDirection(dirx,diry,0.0);
         latDirection.rotateUz(oldDirection);
-
-        G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-        G4double fac = 1.;
-        if(r >  safety) {
-          //  ******* so safety is computed at boundary too ************
-	  G4double newsafety = safetyHelper->ComputeSafety(Position);
-          //G4double newsafety = safety;
-          if(r > newsafety)
-            fac = newsafety/r ;
-        }  
-
-        if(fac > 0.)
-        {
-          // compute new endpoint of the Step
-          G4ThreeVector newPosition = Position+fac*r*latDirection;
-
-	  // definetly not on boundary
-	  if(1. == fac) {
-	    //if(0. < fac) {
-	    safetyHelper->ReLocateWithinVolume(newPosition);
-
-	    
-	  } else {
-            // check safety after displacement
-	    G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	    // displacement to boundary
-            if(postsafety <= 0.) {
-	      safetyHelper->Locate(newPosition, newDirection);
-
-	    // not on the boundary
-            } else { 
-	      safetyHelper->ReLocateWithinVolume(newPosition);
-	    }
-	  }
-          fParticleChange->ProposePosition(newPosition);
-        } 
-     }
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }
   }
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4WaterStopping.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WaterStopping.cc,v 1.11 2008/12/18 13:01:38 gunter Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4WaterStopping.cc,v 1.16 2009/05/15 19:04:21 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 
 //---------------------------------------------------------------------------
@@ -36,4 +36,7 @@
 //
 // Modifications:
+// 29.04.2009 A.Ivantchenko added revised data for G4WATER, provided by 
+//            Prof.P.Sigmund Univ. Southern Denmark in the framework of 
+//            the ESA Technology Research Programme (ESA contract 21435/08/NL/AT)
 //
 //----------------------------------------------------------------------------
@@ -98,5 +101,5 @@
   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[16] = {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};
+  // 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++) {
     Z[i] = zz[i];
@@ -109,51 +112,36 @@
   G4double factor = 1000.*MeV/cm;
 
-G4double G4_WATER_Li [53]={2.626, 2.84, 3.191, 3.461, 3.665, 3.817, 3.927, 4.004, 4.056, 4.102, 3.998, 3.853, 3.702, 3.413, 3.158, 2.934, 2.739, 2.567, 2.415, 2.28, 1.782, 1.465, 1.247, 1.087, 0.8706, 0.7299, 0.631, 0.5575, 0.5005, 0.455, 0.4178, 0.3004, 0.2376, 0.1981, 0.1708, 0.1354, 0.1132, 0.09803, 0.08692, 0.07842, 0.07171, 0.06627, 0.0495, 0.04085, 0.03557, 0.03203, 0.0276, 0.02498, 0.02327, 0.0221, 0.02126, 0.02064, 0.02016, };
-AddData(E,G4_WATER_Li,factor);
-
-G4double G4_WATER_Be [53]={3.272, 3.565, 4.061, 4.463, 4.79, 5.052, 5.258, 5.419, 5.542, 5.803, 5.787, 5.675, 5.529, 5.215, 4.912, 4.634, 4.381, 4.152, 3.944, 3.756, 3.026, 2.533, 2.179, 1.913, 1.542, 1.296, 1.122, 0.9911, 0.8898, 0.8087, 0.7423, 0.5335, 0.4219, 0.3518, 0.3034, 0.2406, 0.2013, 0.1743, 0.1545, 0.1394, 0.1275, 0.1178, 0.08805, 0.07266, 0.06328, 0.05698, 0.0491, 0.04444, 0.04141, 0.03933, 0.03783, 0.03672, 0.03588, };
-AddData(E,G4_WATER_Be,factor);
-
-G4double G4_WATER_B [53]={3.773, 4.142, 4.776, 5.304, 5.749, 6.122, 6.431, 6.684, 6.89, 7.432, 7.551, 7.505, 7.391, 7.091, 6.772, 6.463, 6.172, 5.901, 5.65, 5.418, 4.484, 3.817, 3.322, 2.94, 2.392, 2.02, 1.752, 1.549, 1.391, 1.265, 1.161, 0.8332, 0.6587, 0.5492, 0.4737, 0.3757, 0.3144, 0.2723, 0.2415, 0.2179, 0.1993, 0.1842, 0.1376, 0.1136, 0.09894, 0.08909, 0.07678, 0.0695, 0.06477, 0.06151, 0.05916, 0.05743, 0.05611, };
-AddData(E,G4_WATER_B,factor);
-
-G4double G4_WATER_C [53]={4.154, 4.593, 5.358, 6.009, 6.568, 7.049, 7.46, 7.809, 8.103, 8.968, 9.262, 9.311, 9.25, 8.994, 8.68, 8.358, 8.045, 7.747, 7.465, 7.199, 6.093, 5.269, 4.636, 4.137, 3.403, 2.891, 2.516, 2.23, 2.004, 1.823, 1.673, 1.2, 0.9483, 0.7904, 0.6817, 0.5406, 0.4525, 0.392, 0.3477, 0.3138, 0.287, 0.2653, 0.1983, 0.1637, 0.1426, 0.1284, 0.1107, 0.1002, 0.09335, 0.08865, 0.08528, 0.08278, 0.08088, };
-AddData(E,G4_WATER_C,factor);
-
-G4double G4_WATER_N [53]={4.49, 4.984, 5.86, 6.616, 7.276, 7.854, 8.36, 8.799, 9.179, 10.39, 10.89, 11.07, 11.08, 10.9, 10.61, 10.3, 9.974, 9.66, 9.357, 9.068, 7.823, 6.859, 6.097, 5.484, 4.56, 3.9, 3.408, 3.029, 2.727, 2.482, 2.28, 1.636, 1.291, 1.076, 0.9274, 0.7354, 0.6156, 0.5333, 0.4731, 0.427, 0.3906, 0.3611, 0.27, 0.2229, 0.1942, 0.1749, 0.1507, 0.1365, 0.1272, 0.1208, 0.1162, 0.1128, 0.1102, };
-AddData(E,G4_WATER_N,factor);
-
-G4double G4_WATER_O [53]={4.778, 5.321, 6.298, 7.152, 7.907, 8.578, 9.173, 9.7, 10.16, 11.73, 12.46, 12.78, 12.89, 12.81, 12.57, 12.27, 11.95, 11.63, 11.32, 11.01, 9.659, 8.571, 7.691, 6.967, 5.854, 5.042, 4.427, 3.945, 3.56, 3.245, 2.983, 2.142, 1.689, 1.406, 1.212, 0.9602, 0.8037, 0.6963, 0.6178, 0.5577, 0.5102, 0.4716, 0.3527, 0.2913, 0.2538, 0.2285, 0.197, 0.1784, 0.1663, 0.1579, 0.1519, 0.1475, 0.1441, };
-AddData(E,G4_WATER_O,factor);
-
-G4double G4_WATER_F [53]={4.992, 5.575, 6.637, 7.578, 8.418, 9.171, 9.847, 10.45, 11, 12.9, 13.88, 14.35, 14.56, 14.59, 14.4, 14.13, 13.83, 13.51, 13.19, 12.87, 11.44, 10.26, 9.279, 8.463, 7.187, 6.237, 5.506, 4.928, 4.461, 4.076, 3.753, 2.707, 2.137, 1.779, 1.533, 1.215, 1.017, 0.8809, 0.7816, 0.7056, 0.6456, 0.5969, 0.4466, 0.3688, 0.3213, 0.2894, 0.2496, 0.2259, 0.2106, 0.2, 0.1924, 0.1868, 0.1825, };
-AddData(E,G4_WATER_F,factor);
-
-G4double G4_WATER_Ne [53]={5.182, 5.797, 6.931, 7.948, 8.865, 9.693, 10.44, 11.12, 11.74, 13.98, 15.21, 15.85, 16.17, 16.33, 16.21, 15.98, 15.69, 15.38, 15.06, 14.74, 13.24, 11.98, 10.91, 10.01, 8.584, 7.503, 6.66, 5.986, 5.436, 4.979, 4.595, 3.332, 2.635, 2.195, 1.892, 1.499, 1.255, 1.087, 0.9646, 0.8709, 0.7969, 0.7368, 0.5514, 0.4555, 0.3969, 0.3576, 0.3083, 0.2792, 0.2602, 0.2472, 0.2378, 0.2308, 0.2255, };
-AddData(E,G4_WATER_Ne,factor);
-
-G4double G4_WATER_Na [53]={5.352, 5.998, 7.203, 8.3, 9.298, 10.21, 11.04, 11.81, 12.5, 15.13, 16.68, 17.56, 18.04, 18.43, 18.44, 18.29, 18.05, 17.78, 17.48, 17.18, 15.67, 14.32, 13.15, 12.14, 10.5, 9.226, 8.218, 7.401, 6.728, 6.166, 5.69, 4.112, 3.237, 2.686, 2.307, 1.821, 1.521, 1.317, 1.168, 1.054, 0.9644, 0.8917, 0.6674, 0.5514, 0.4806, 0.4329, 0.3734, 0.3381, 0.3152, 0.2993, 0.288, 0.2796, 0.2732, };
-AddData(E,G4_WATER_Na,factor);
-
-G4double G4_WATER_Mg [53]={5.542, 6.193, 7.42, 8.551, 9.59, 10.54, 11.42, 12.23, 12.98, 15.85, 17.62, 18.66, 19.26, 19.76, 19.83, 19.7, 19.47, 19.2, 18.89, 18.58, 17.02, 15.62, 14.41, 13.36, 11.64, 10.3, 9.233, 8.362, 7.64, 7.033, 6.516, 4.777, 3.792, 3.162, 2.725, 2.159, 1.806, 1.565, 1.388, 1.254, 1.147, 1.061, 0.7944, 0.6565, 0.5722, 0.5156, 0.4447, 0.4027, 0.3754, 0.3566, 0.343, 0.333, 0.3254, };
-AddData(E,G4_WATER_Mg,factor);
-
-G4double G4_WATER_Al [53]={5.724, 6.39, 7.649, 8.82, 9.905, 10.91, 11.84, 12.71, 13.51, 16.66, 18.69, 19.93, 20.68, 21.38, 21.56, 21.5, 21.32, 21.07, 20.79, 20.48, 18.91, 17.47, 16.19, 15.08, 13.24, 11.78, 10.61, 9.641, 8.835, 8.153, 7.569, 5.583, 4.444, 3.71, 3.199, 2.534, 2.12, 1.836, 1.629, 1.471, 1.346, 1.245, 0.9325, 0.7707, 0.6719, 0.6054, 0.5223, 0.473, 0.441, 0.4189, 0.403, 0.3912, 0.3823, };
-AddData(E,G4_WATER_Al,factor);
-
-G4double G4_WATER_Si [53]={5.905, 6.583, 7.868, 9.073, 10.2, 11.25, 12.23, 13.14, 13.99, 17.4, 19.66, 21.1, 22.01, 22.91, 23.21, 23.22, 23.09, 22.87, 22.61, 22.32, 20.76, 19.28, 17.95, 16.78, 14.83, 13.26, 11.99, 10.94, 10.06, 9.304, 8.656, 6.43, 5.135, 4.294, 3.705, 2.938, 2.458, 2.129, 1.889, 1.706, 1.561, 1.444, 1.082, 0.8942, 0.7796, 0.7026, 0.6061, 0.549, 0.5119, 0.4862, 0.4678, 0.4542, 0.4438, };
-AddData(E,G4_WATER_Si,factor);
-
-G4double G4_WATER_P [53]={6.12, 6.81, 8.118, 9.352, 10.51, 11.61, 12.63, 13.58, 14.48, 18.13, 20.63, 22.28, 23.34, 24.47, 24.91, 25.02, 24.95, 24.78, 24.55, 24.28, 22.76, 21.26, 19.89, 18.67, 16.59, 14.92, 13.54, 12.39, 11.42, 10.59, 9.867, 7.367, 5.896, 4.935, 4.259, 3.376, 2.824, 2.445, 2.169, 1.959, 1.792, 1.657, 1.242, 1.027, 0.8954, 0.807, 0.6963, 0.6308, 0.5881, 0.5587, 0.5375, 0.5219, 0.51, };
-AddData(E,G4_WATER_P,factor);
-
-G4double G4_WATER_S [53]={6.294, 7, 8.338, 9.604, 10.8, 11.94, 13, 14, 14.94, 18.82, 21.55, 23.41, 24.65, 26.01, 26.6, 26.81, 26.81, 26.69, 26.5, 26.26, 24.79, 23.28, 21.88, 20.61, 18.43, 16.64, 15.16, 13.92, 12.86, 11.95, 11.15, 8.371, 6.715, 5.624, 4.856, 3.847, 3.217, 2.785, 2.47, 2.229, 2.04, 1.886, 1.413, 1.169, 1.019, 0.9187, 0.7929, 0.7183, 0.6697, 0.6362, 0.6122, 0.5944, 0.5808, };
-AddData(E,G4_WATER_S,factor);
-
-G4double G4_WATER_Cl [53]={6.522, 7.237, 8.59, 9.875, 11.1, 12.26, 13.36, 14.39, 15.37, 19.45, 22.4, 24.45, 25.86, 27.46, 28.19, 28.5, 28.57, 28.5, 28.34, 28.13, 26.72, 25.21, 23.78, 22.47, 20.2, 18.32, 16.75, 15.42, 14.28, 13.3, 12.44, 9.392, 7.557, 6.34, 5.479, 4.344, 3.633, 3.145, 2.789, 2.517, 2.303, 2.13, 1.596, 1.32, 1.151, 1.038, 0.8957, 0.8115, 0.7567, 0.7189, 0.6917, 0.6717, 0.6563, };
-AddData(E,G4_WATER_Cl,factor);
-
-G4double G4_WATER_Ar [53]={6.642, 7.369, 8.739, 10.04, 11.28, 12.47, 13.59, 14.66, 15.66, 19.93, 23.08, 25.32, 26.89, 28.72, 29.6, 30.01, 30.15, 30.13, 30, 29.82, 28.46, 26.94, 25.49, 24.15, 21.81, 19.86, 18.22, 16.82, 15.61, 14.57, 13.65, 10.39, 8.394, 7.063, 6.114, 4.859, 4.067, 3.522, 3.125, 2.821, 2.581, 2.387, 1.789, 1.48, 1.291, 1.164, 1.005, 0.9105, 0.8491, 0.8067, 0.7763, 0.7537, 0.7365, };
-AddData(E,G4_WATER_Ar,factor);
+  G4double G4_WATER_Li[53]={2.3193,2.5198,2.8539,3.1164,3.3203,3.4756,3.5914,3.6755,3.7347,3.8125,3.7349,3.6134,3.4818,3.2258,2.9949,2.7909,2.611,2.4517,2.3103,2.1841,1.7151,1.4139,1.2053,1.0525,0.84417,0.70862,0.61317,0.54214,0.48708,0.44305,0.40697,0.29312,0.23208,0.19364,0.16706,0.13252,0.11092,0.09608,0.08522,0.076915,0.07035,0.065026,0.048615,0.040137,0.034964,0.03149,0.027148,0.024579,0.022911,0.021761,0.020937,0.020327,0.019862};
+  AddData(E,G4_WATER_Li,factor);
+  G4double G4_WATER_Be[53]={2.872,3.1439,3.6102,3.9967,4.3169,4.5788,4.7897,4.9568,5.0872,5.387,5.4028,5.3185,5.1971,4.9243,4.6549,4.4036,4.1732,3.9629,3.771,3.5957,2.9117,2.4439,2.1062,1.8518,1.4956,1.2587,1.0901,0.96393,0.86589,0.78742,0.72313,0.52053,0.41214,0.34394,0.2968,0.2355,0.19717,0.17081,0.15152,0.13676,0.1251,0.11564,0.086471,0.071399,0.062202,0.056023,0.048303,0.043734,0.040767,0.038723,0.037256,0.036172,0.035345};
+  AddData(E,G4_WATER_Be,factor);
+  G4double G4_WATER_B[53]={3.2922,3.6315,4.2226,4.7242,5.1543,5.5219,5.8323,6.0911,6.3043,6.8888,7.0451,7.0309,6.9445,6.6925,6.4129,6.1372,5.8747,5.6283,5.3983,5.1841,4.3121,3.6826,3.2109,2.8459,2.3203,1.9619,1.7028,1.5072,1.3543,1.2315,1.1307,0.81305,0.64344,0.53693,0.46337,0.36777,0.30797,0.26684,0.23674,0.21371,0.1955,0.18072,0.13517,0.11163,0.097256,0.087601,0.075535,0.068395,0.063757,0.060561,0.058268,0.056574,0.05528};
+  AddData(E,G4_WATER_B,factor);
+  G4double G4_WATER_C[53]={3.6037,4.0035,4.7125,5.3248,5.8601,6.3287,6.7363,7.0875,7.3872,8.2986,8.635,8.7189,8.6879,8.485,8.2162,7.9331,7.6534,7.3839,7.1272,6.884,5.8573,5.0814,4.4808,4.0044,3.3005,2.808,2.4458,2.1691,1.9511,1.7751,1.6302,1.1714,0.9263,0.77269,0.66676,0.52925,0.44328,0.38415,0.34086,0.30773,0.28154,0.26028,0.19473,0.16083,0.14014,0.12624,0.10886,0.098575,0.091894,0.08729,0.083986,0.081545,0.079682};
+  AddData(E,G4_WATER_C,factor);
+  G4double G4_WATER_N[53]={3.8821,4.3278,5.1314,5.8371,6.4632,7.021,7.5168,7.9545,8.3378,9.5943,10.145,10.356,10.402,10.278,10.041,9.7677,9.4845,9.2035,8.9301,8.6668,7.5173,6.6126,5.8919,5.308,4.4226,3.7883,3.3138,2.9467,2.655,2.4179,2.2217,1.5965,1.2614,1.0516,0.90715,0.71995,0.60305,0.52268,0.46384,0.41882,0.3832,0.35431,0.26515,0.21903,0.19087,0.17194,0.14829,0.13429,0.12519,0.11892,0.11442,0.1111,0.10856};
+  AddData(E,G4_WATER_N,factor);
+  G4double G4_WATER_O[53]={4.1215,4.6063,5.4945,6.2868,6.9978,7.6391,8.2175,8.7372,9.201,10.808,11.596,11.955,12.096,12.077,11.89,11.639,11.364,11.081,10.799,10.523,9.2787,8.2615,7.4307,6.7431,5.6787,4.8981,4.3045,3.8393,3.4663,3.161,2.9069,2.0903,1.6501,1.3745,1.1851,0.94004,0.78733,0.68244,0.60568,0.54694,0.50048,0.46278,0.34643,0.28622,0.24945,0.22474,0.19384,0.17555,0.16367,0.15547,0.14959,0.14525,0.14194};
+  AddData(E,G4_WATER_O,factor);
+  G4double G4_WATER_F[53]={4.2951,4.8107,5.7678,6.6342,7.4196,8.1343,8.7858,9.3786,9.9152,11.857,12.89,13.408,13.652,13.749,13.62,13.398,13.136,12.857,12.573,12.291,10.982,9.88,8.9601,8.1871,6.9687,6.0574,5.3535,4.7954,4.3434,3.9704,3.658,2.642,2.0878,1.7393,1.4994,1.1892,0.99601,0.86336,0.76632,0.69207,0.63334,0.58568,0.43857,0.36242,0.3159,0.28462,0.24552,0.22237,0.20733,0.19696,0.18951,0.18401,0.17981};
+  AddData(E,G4_WATER_F,factor);
+  G4double G4_WATER_Ne[53]={4.4513,4.991,6.004,6.9338,7.7852,8.5662,9.284,9.9435,10.547,12.813,14.099,14.791,15.151,15.382,15.321,15.136,14.895,14.626,14.345,14.061,12.705,11.53,10.532,9.6823,8.3208,7.2846,6.4735,5.8234,5.2919,4.8502,4.4779,3.252,2.5745,2.146,1.8503,1.4675,1.2291,1.0654,0.94575,0.85419,0.78178,0.72301,0.54158,0.44763,0.39022,0.35162,0.30335,0.27477,0.25619,0.24339,0.23419,0.2274,0.22222};
+  AddData(E,G4_WATER_Ne,factor);
+  G4double G4_WATER_Na[53]={4.5914,5.1553,6.2244,7.2203,8.1435,8.9982,9.7906,10.526,11.206,13.848,15.453,16.384,16.916,17.369,17.442,17.344,17.16,16.93,16.675,16.407,15.045,13.799,12.706,11.753,10.187,8.9672,7.9956,7.2072,6.5562,6.0112,5.5493,4.0154,3.1635,2.6262,2.2573,1.7832,1.4904,1.2907,1.1451,1.0339,0.94615,0.87498,0.65548,0.54186,0.47243,0.42574,0.36734,0.33275,0.31027,0.29478,0.28364,0.27542,0.26915};
+  AddData(E,G4_WATER_Na,factor);
+  G4double G4_WATER_Mg[53]={4.7537,5.3178,6.3962,7.4137,8.3663,9.2554,10.085,10.859,11.581,14.455,16.279,17.373,18.018,18.598,18.727,18.654,18.479,18.25,17.99,17.716,16.313,15.026,13.895,12.907,11.277,9.9981,8.9727,8.1344,7.4376,6.8507,6.35,4.6625,3.7049,3.0916,2.666,2.1135,1.7695,1.5336,1.3613,1.2296,1.1255,1.041,0.7802,0.64511,0.56252,0.50698,0.43749,0.39634,0.36958,0.35114,0.33789,0.3281,0.32063};
+  AddData(E,G4_WATER_Mg,factor);
+  G4double G4_WATER_Al[53]={4.911,5.4856,6.5852,7.6302,8.6193,9.551,10.426,11.248,12.018,15.152,17.23,18.531,19.33,20.11,20.354,20.352,20.224,20.024,19.785,19.521,18.12,16.795,15.611,14.565,12.819,11.431,10.305,9.3772,8.6003,7.9413,7.3761,5.45,4.342,3.6272,3.129,2.4808,2.0767,1.7997,1.5975,1.4429,1.3207,1.2216,0.91581,0.75739,0.66053,0.59537,0.51383,0.46554,0.43413,0.41248,0.39693,0.38544,0.37667};
+  AddData(E,G4_WATER_Al,factor);
+  G4double G4_WATER_Si[53]={5.0697,5.6529,6.77,7.8376,8.8564,9.8229,10.736,11.597,12.409,15.773,18.09,19.594,20.549,21.535,21.9,21.975,21.896,21.73,21.513,21.265,19.878,18.525,17.298,16.202,14.354,12.867,11.651,10.64,9.7876,9.0607,8.4341,6.2756,5.0169,4.1981,3.6247,2.8758,2.4078,2.0868,1.8523,1.6731,1.5315,1.4167,1.0623,0.8787,0.76644,0.69091,0.59637,0.54036,0.50394,0.47883,0.46078,0.44746,0.43728};
+  AddData(E,G4_WATER_Si,factor);
+  G4double G4_WATER_P[53]={5.2616,5.8538,6.9867,8.074,9.1192,10.117,11.065,11.964,12.815,16.396,18.945,20.657,21.779,22.993,23.502,23.672,23.659,23.54,23.356,23.132,21.792,20.425,19.162,18.02,16.066,14.472,13.155,12.051,11.114,10.311,9.6149,7.1917,5.7614,4.8249,4.1667,3.305,2.7661,2.3966,2.1269,1.9215,1.7603,1.6263,1.2197,1.0091,0.88027,0.79361,0.68512,0.62083,0.57902,0.55019,0.52947,0.51417,0.50248};
+  AddData(E,G4_WATER_P,factor);
+  G4double G4_WATER_S[53]={5.4129,6.0193,7.1761,8.2871,9.36,10.39,11.373,12.308,13.196,16.986,19.762,21.683,22.976,24.431,25.091,25.363,25.421,25.354,25.208,25.013,23.734,22.366,21.074,19.891,17.84,16.146,14.731,13.536,12.516,11.636,10.869,8.1727,6.5617,5.4997,4.7505,3.7671,3.1514,2.7293,2.4215,2.1866,2.0012,1.8509,1.3881,1.1485,1.002,0.90348,0.78008,0.70695,0.65938,0.62657,0.60299,0.58558,0.57229};
+  AddData(E,G4_WATER_S,factor);
+  G4double G4_WATER_Cl[53]={5.6171,6.2307,7.3984,8.5209,9.6097,10.661,11.669,12.632,13.551,17.518,20.497,22.615,24.076,25.769,26.58,26.953,27.082,27.066,26.958,26.791,25.579,24.214,22.901,21.685,19.553,17.771,16.271,14.996,13.9,12.949,12.118,9.1688,7.385,6.2003,5.3604,4.2535,3.5588,3.0821,2.7343,2.4689,2.2595,2.0898,1.5673,1.297,1.1318,1.0205,0.88128,0.79874,0.74504,0.708,0.68138,0.66172,0.64671};
+  AddData(E,G4_WATER_Cl,factor);
+  G4double G4_WATER_Ar[53]={5.7158,6.3394,7.5204,8.6525,9.7528,10.82,11.849,12.836,13.78,17.904,21.07,23.375,24.999,26.928,27.889,28.361,28.559,28.592,28.521,28.381,27.228,25.869,24.541,23.299,21.103,19.255,17.686,16.346,15.19,14.183,13.3,10.137,8.2021,6.9062,5.9819,4.757,3.9841,3.4523,3.0636,2.7668,2.5324,2.3425,1.7573,1.4546,1.2694,1.1448,0.98872,0.8962,0.83601,0.79448,0.76464,0.74259,0.72575};
+  AddData(E,G4_WATER_Ar,factor);
 
   if(corr) {
Index: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4WentzelVIModel.cc,v 1.16 2008/11/19 11:47:50 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4WentzelVIModel.cc,v 1.32 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -52,5 +52,4 @@
 //
 
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -60,6 +59,4 @@
 #include "G4LossTableManager.hh"
 #include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
 #include "G4PhysicsTableHelper.hh"
 #include "G4ElementVector.hh"
@@ -71,4 +68,7 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
+G4double G4WentzelVIModel::FormFactor[]    = {0.0};
 
 using namespace std;
@@ -96,8 +96,7 @@
   thePositron = G4Positron::Positron();
   theProton   = G4Proton::Proton();
-  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
+  lowEnergyLimit = 0.1*keV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff  = twopi*p0*p0;
-  constn = 6.937e-6/(MeV*MeV);
   tkin = targetZ = mom2 = DBL_MIN;
   ecut = etag = DBL_MAX;
@@ -106,7 +105,20 @@
   xsecn.resize(nelments);
   prob.resize(nelments);
-  for(size_t j=0; j<100; j++) {
-    FF[j]    = 0.0;
-  } 
+
+  // Thomas-Fermi screening radii
+  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
+
+  if(0.0 == ScreenRSquare[0]) {
+    G4double a0 = electron_mass_c2/0.88534; 
+    G4double constn = 6.937e-6/(MeV*MeV);
+
+    ScreenRSquare[0] = alpha2*a0*a0;
+    for(G4int j=1; j<100; j++) {
+      G4double x = a0*fNistManager->GetZ13(j);
+      ScreenRSquare[j] = alpha2*x*x;
+      x = fNistManager->GetA27(j); 
+      FormFactor[j] = constn*x*x;
+    } 
+  }
 }
 
@@ -132,13 +144,6 @@
   if(!isInitialized) {
     isInitialized = true;
-
-    if (pParticleChange)
-      fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-    else
-      fParticleChange = new G4ParticleChangeForMSC();
-
-    safetyHelper = G4TransportationManager::GetTransportationManager()
-      ->GetSafetyHelper();
-    safetyHelper->InitialiseHelper();
+    fParticleChange = GetParticleChangeForMSC();
+    InitialiseSafetyHelper();
   }
 }
@@ -156,5 +161,5 @@
   SetupKinematic(ekin, cutEnergy);
   SetupTarget(Z, ekin);
-  G4double xsec = ComputeTransportXSectionPerVolume();
+  G4double xsec = ComputeTransportXSectionPerAtom();
   /*  
   G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2 
@@ -167,5 +172,5 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume()
+G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
 {
   G4double xSection = 0.0;
@@ -189,26 +194,27 @@
       y = 0.0;
     }
-    xSection += y/targetZ;
-  }
-  /*
-  G4cout << "G4WentzelVIModel:XS per A " << " Z= " << Z << " e(MeV)= " << kinEnergy/MeV 
+    xSection = y;
+  }
+  /*  
+  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ 
+	 << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
 	 << " cut(MeV)= " << ecut/MeV  
   	 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ 
-	 << " zmaxN= " << (1.0 - cosTetMsxNuc)/screenZ << G4endl;
+	 << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ 
+         << " costm= " << cosTetMaxNuc2 << G4endl;
   */
-
   // scattering off nucleus
   if(cosTetMaxNuc2 < 1.0) {
     x  = 1.0 - cosTetMaxNuc2;
     x1 = screenZ*formfactA;
-    x2 = 1.0/(1.0 - x1); 
+    x2 = 1.0 - x1; 
     x3 = x/screenZ;
     x4 = formfactA*x;
     // low-energy limit
     if(x3 < numlimit && x1 < numlimit) {
-      y = 0.5*x3*x3*x2*x2*x2*(1.0 - 1.333333*x3 + 1.5*x3*x3 - 1.5*x1
-      			      + 3.0*x1*x1 + 2.666666*x3*x1);
+      y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1
+		     + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2);
       // high energy limit
-    } else if(1.0 < x1) {
+    } else if(x2 <= 0.0) {
       x4 = x1*(1.0 + x3);
       y  = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
@@ -216,6 +222,8 @@
     } else {
       y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4)) 
-	   - x3/(1. + x3) - x4/(1. + x4))*x2*x2; 
-    }
+	   - x3/(1. + x3) - x4/(1. + x4))/(x2*x2); 
+    }
+    //G4cout << "y= " << y << " x1= " <<x1<<"  x2= " <<x2
+    //	   <<"  x3= "<<x3<<"  x4= " << x4<<G4endl;
     if(y < 0.0) {
       nwarnings++;
@@ -230,9 +238,15 @@
       y = 0.0;
     }
-    xSection += y; 
-  }
-  xSection *= (coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2); 
-  //  G4cout << "   XStotal= " << xSection/barn << " screenZ= " << screenZ 
-  //	 << " formF= " << formfactA << " for " << p->GetParticleName() << G4endl;
+    xSection += y*targetZ; 
+  }
+  xSection *= kinFactor;
+  /*
+  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn 
+	 << " screenZ= " << screenZ << " formF= " << formfactA 
+	 << " for " << particle->GetParticleName() 
+  	 << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
+	 << " x= " << x 
+	 << G4endl;
+  */
   return xSection; 
 }
@@ -277,5 +291,5 @@
   // i.e. when it is needed for optimization purposes
   if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 
-    presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
+    presafety = ComputeSafety(sp->GetPosition(), tlimit); 
   /*
   G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= " 
@@ -291,5 +305,6 @@
     G4double rlimit = facrange*lambda0;
     G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
-    if(rcut > rlimit) rlimit = std::pow(2.0*rcut*rcut*lambda0,0.33333333);
+    if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
+    rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
     if(rlimit < tlimit) tlimit = rlimit;
   }
@@ -404,5 +419,5 @@
   } else {
 
-    // define threshold angle as 2 sigma of central value
+    // define threshold angle between single and multiple scattering 
     cosThetaMin = 1.0 - 3.0*x1;
 
@@ -570,35 +585,6 @@
 
     if(r > tlimitminfix) {
-      G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-      G4double fac= 1.;
-      if(r >= safety) {
-	//  ******* so safety is computed at boundary too ************
-	G4double newsafety = 
-	  safetyHelper->ComputeSafety(Position) - tlimitminfix;
-        if(newsafety <= 0.0) fac = 0.0;
-	else if(r > newsafety) fac = newsafety/r ;
-        //G4cout << "NewSafety= " << newsafety << " fac= " << fac 
-	// << " r= " << r << " sint= " << sint << " pos " << Position << G4endl;
-      }  
-
-      if(fac > 0.) {
-	// compute new endpoint of the Step
-	G4ThreeVector newPosition = Position + fac*pos;
-
-	// check safety after displacement
-	G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	// displacement to boundary
-	if(postsafety <= 0.0) {
-	  safetyHelper->Locate(newPosition, newDirection);
-
-	  // not on the boundary
-	} else { 
-	  safetyHelper->ReLocateWithinVolume(newPosition);
-	  // if(fac < 1.0) G4cout << "NewPosition " << newPosition << G4endl;
-	}
-     
-	fParticleChange->ProposePosition(newPosition);
-      } 
+      pos /= r;
+      ComputeDisplacement(fParticleChange, pos, r, safety);
     }
   }
@@ -624,6 +610,4 @@
   G4double xs = 0.0;
 
-  G4double fac = coeff*chargeSquare*invbeta2/mom2;
-
   for (G4int i=0; i<nelm; i++) {
     SetupTarget((*theElementVector)[i]->GetZ(), tkin);
@@ -635,5 +619,5 @@
     cosTetMaxNuc2  = std::max(cosnm,cosThetaMin); 
     cosTetMaxElec2 = std::max(cosem,cosThetaMin); 
-    xtsec += ComputeTransportXSectionPerVolume()*density;
+    xtsec += ComputeTransportXSectionPerAtom()*density;
     // return limit back
     cosTetMaxElec2 = cosem;
@@ -643,5 +627,5 @@
     G4double nsec = 0.0;
     G4double x1 = 1.0 - cosThetaMin + screenZ;
-    G4double f  = fac*targetZ*density; 
+    G4double f  = kinFactor*density; 
 
     // scattering off electrons
@@ -656,5 +640,6 @@
       G4double s  = screenZ*formfactA;
       G4double z1 = 1.0 - cosnm + screenZ;
-      G4double d  = (1.0 - s)/formfactA;
+      G4double s1 = 1.0 - s;
+      G4double d  = s1/formfactA;
 
       // check numerical limit
@@ -667,6 +652,5 @@
 	G4double x2 = x1 + d;
 	G4double z2 = z1 + d;
-	nsec = (1.0 + 2.0*s)*((cosThetaMin - cosnm)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
-			   2.0*log(z1*x2/(z2*x1))/d);
+	nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
       }
       nsec *= f*targetZ;
@@ -801,5 +785,6 @@
       G4double mom22 = t1*(t1 + 2.0*mass);
       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
-      if(ctm < 1.0) cosTetMaxElec = ctm;
+      if(ctm <  1.0) cosTetMaxElec = ctm;
+      if(ctm < -1.0) cosTetMaxElec = -1.0;
     }
   }
@@ -808,11 +793,2 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4WentzelVIModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
-                                     const G4MaterialCutsCouple*,
-				     const G4DynamicParticle*,
-				     G4double,
-				     G4double)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlung.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlung.cc,v 1.55 2008/11/14 19:23:07 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlung.cc,v 1.56 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -104,4 +104,11 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4bool G4eBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4eBremsstrahlung::InitialiseEnergyLossProcess(
                                                 const G4ParticleDefinition* p,
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlungModel.cc,v 1.43 2008/11/13 19:28:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlungModel.cc,v 1.44 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -160,10 +160,5 @@
   }
   if(isInitialised) return;
-
-  if(pParticleChange) {
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-  } else {
-    fParticleChange = new G4ParticleChangeForLoss();
-  }
+  fParticleChange = GetParticleChangeForLoss();
   isInitialised = true;
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eBremsstrahlungRelModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eBremsstrahlungRelModel.cc,v 1.12 2008/11/13 23:28:27 schaelic Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eBremsstrahlungRelModel.cc,v 1.14 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -178,10 +178,5 @@
 
   if(isInitialised) return;
-
-  if(pParticleChange) {
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-  } else {
-    fParticleChange = new G4ParticleChangeForLoss();
-  }
+  fParticleChange = GetParticleChangeForLoss();
   isInitialised = true;
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eCoulombScatteringModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eCoulombScatteringModel.cc,v 1.59 2008/10/22 18:39:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eCoulombScatteringModel.cc,v 1.69 2009/05/10 16:09:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -66,4 +66,7 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4double G4eCoulombScatteringModel::ScreenRSquare[] = {0.0};
+G4double G4eCoulombScatteringModel::FormFactor[]    = {0.0};
+
 using namespace std;
 
@@ -84,17 +87,29 @@
   currentMaterial = 0; 
   currentElement  = 0;
-  a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
+  lowEnergyLimit = keV;
   G4double p0 = electron_mass_c2*classic_electr_radius;
   coeff  = twopi*p0*p0;
-  constn = 6.937e-6/(MeV*MeV);
   tkin = targetZ = mom2 = DBL_MIN;
   elecXSection = nucXSection = 0.0;
-  recoilThreshold = DBL_MAX;
+  recoilThreshold = 100.*keV;
   ecut = DBL_MAX;
   particle = 0;
   currentCouple = 0;
-  for(size_t j=0; j<100; j++) {
-    FF[j] = 0.0;
-  } 
+
+  // Thomas-Fermi screening radii
+  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
+
+  if(0.0 == ScreenRSquare[0]) {
+    G4double a0 = electron_mass_c2/0.88534; 
+    G4double constn = 6.937e-6/(MeV*MeV);
+
+    ScreenRSquare[0] = alpha2*a0*a0;
+    for(G4int j=1; j<100; j++) {
+      G4double x = a0*fNistManager->GetZ13(j);
+      ScreenRSquare[j] = alpha2*x*x;
+      x = fNistManager->GetA27(j); 
+      FormFactor[j] = constn*x*x;
+    } 
+  }
 }
 
@@ -121,10 +136,5 @@
   if(!isInitialised) {
     isInitialised = true;
-
-    if(pParticleChange)
-      fParticleChange = 
-	reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-    else
-      fParticleChange = new G4ParticleChangeForGamma();
+    fParticleChange = GetParticleChangeForGamma();
   }
   if(mass < GeV && particle->GetParticleType() != "nucleus") {
@@ -157,5 +167,6 @@
       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
       //G4cout << "ctm= " << ctm << G4endl;
-      if(ctm < 1.0) cosTetMaxElec = ctm;
+      if(ctm <  1.0) cosTetMaxElec = ctm;
+      if(ctm < -1.0) cosTetMaxElec = -1.0;
     }
   }
@@ -196,6 +207,5 @@
 {
   // This method needs initialisation before be called
-
-  G4double fac = coeff*targetZ*chargeSquare*invbeta2/mom2; 
+  G4double fac = coeff*targetZ*chargeSquare*kinFactor;
   elecXSection = 0.0;
   nucXSection  = 0.0;
@@ -216,5 +226,6 @@
     G4double s  = screenZ*formfactA;
     G4double z1 = 1.0 - cosTetMaxNuc2 + screenZ;
-    G4double d  = (1.0 - s)/formfactA;
+    G4double s1 = 1.0 - s;
+    G4double d  = s1/formfactA;
     //G4cout <<"x1= "<<x1<<" z1= " <<z1<<" s= "<<s << " d= " <<d <<G4endl;
     if(d < 0.2*x1) {
@@ -226,10 +237,8 @@
       G4double x2 = x1 + d;
       G4double z2 = z1 + d;
-      x = (1.0 + 2.0*s)*((cosTetMinNuc - cosTetMaxNuc2)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
-			 2.0*log(z1*x2/(z2*x1))/d);
+      x = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
     }
     nucXSection += fac*targetZ*x;
   }
-
   //G4cout<<" cross(bn)= "<<nucXSection/barn<<" xsElec(bn)= "<<elecXSection/barn
   //	<< " Asc= " << screenZ << G4endl; 
@@ -283,10 +292,6 @@
     G4int ia = SelectIsotopeNumber(currentElement);
     G4double Trec = z1*mom2/(amu_c2*G4double(ia));
-    G4double th = 
-      std::min(recoilThreshold,
-	       targetZ*currentElement->GetIonisation()->GetMeanExcitationEnergy());
-
-    if(Trec > th) {
-      G4int iz = G4int(targetZ);
+
+    if(Trec > recoilThreshold) {
       G4ParticleDefinition* ion = theParticleTable->FindIon(iz, ia, 0, iz);
       Trec = z1*mom2/ion->GetPDGMass();
Index: trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eIonisation.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eIonisation.cc,v 1.56 2008/10/20 08:56:41 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eIonisation.cc,v 1.57 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -103,4 +103,22 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4double G4eIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
+					 const G4Material*,
+					 G4double cut)
+{
+  G4double x = cut;
+  if(isElectron) x += cut;
+  return x;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4eIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Electron::Electron() || &p == G4Positron::Positron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4eIonisation::InitialiseEnergyLossProcess(
 		    const G4ParticleDefinition* part,
Index: trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eeToTwoGammaModel.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eeToTwoGammaModel.cc,v 1.14 2007/05/23 08:47:35 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eeToTwoGammaModel.cc,v 1.15 2009/04/09 18:41:18 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -102,11 +102,5 @@
 {
   if(isInitialised) return;
-
-  if(pParticleChange)
-    fParticleChange =
-      reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForGamma();
-
+  fParticleChange = GetParticleChangeForGamma();
   isInitialised = true;
 }
Index: trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4eplusAnnihilation.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4eplusAnnihilation.cc,v 1.29 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4eplusAnnihilation.cc,v 1.30 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -79,4 +79,20 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4bool G4eplusAnnihilation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (&p == G4Positron::Positron());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4eplusAnnihilation::AtRestGetPhysicalInteractionLength(
+                              const G4Track&, G4ForceCondition* condition)
+{
+  *condition = NotForced;
+  return 0.0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4eplusAnnihilation::InitialiseProcess(const G4ParticleDefinition*)
 {
Index: trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4hIonisation.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4hIonisation.cc,v 1.81 2008/10/22 16:02:20 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4hIonisation.cc,v 1.82 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -122,4 +122,23 @@
 {}
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4hIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV &&
+	 !p.IsShortLived());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4hIonisation::MinPrimaryEnergy(const G4ParticleDefinition*,
+					 const G4Material*,
+					 G4double cut)
+{
+  G4double x = 0.5*cut/electron_mass_c2;
+  G4double g = x*ratio + std::sqrt((1. + x)*(1. + x*ratio*ratio));
+  return mass*(g - 1.0);
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....  
 
Index: trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc	(revision 1007)
+++ trunk/source/processes/electromagnetic/standard/src/G4ionIonisation.cc	(revision 1055)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4ionIonisation.cc,v 1.65 2008/10/15 17:53:44 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-02 $
+// $Id: G4ionIonisation.cc,v 1.66 2009/02/20 12:06:37 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
 //
 // -------------------------------------------------------------------
@@ -104,4 +104,22 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+G4bool G4ionIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived() &&
+          p.GetParticleType() == "nucleus");
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ionIonisation::MinPrimaryEnergy(const G4ParticleDefinition* p, 
+					   const G4Material*, 
+					   G4double cut)
+{
+  return 
+    p->GetPDGMass()*(std::sqrt(1. + 0.5*cut/CLHEP::electron_mass_c2) - 1.0);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4ionIonisation::InitialiseEnergyLossProcess(
 		      const G4ParticleDefinition* part,
