Index: trunk/source/processes/electromagnetic/standard/include/G4ICRU49NuclearStoppingModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4ICRU49NuclearStoppingModel.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4ICRU49NuclearStoppingModel.hh	(revision 1197)
@@ -0,0 +1,119 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4ICRU49NuclearStoppingModel.hh,v 1.2 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+//
+// GEANT4 Class header file
+//
+//
+// File name:     G4ICRU49NuclearStoppingModel
+//
+// Author:        V.Ivanchenko 
+//
+// Creation date: 20.07.2009 
+//
+// Modifications:
+//
+//
+// Class Description:
+//
+// Implementation of the ICRU'49 model of nuclear stopping 
+
+// -------------------------------------------------------------------
+//
+
+#ifndef G4ICRU49NuclearStoppingModel_h
+#define G4ICRU49NuclearStoppingModel_h 1
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4VEmModel.hh"
+
+class G4ParticleChangeForLoss;
+class G4Pow;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4ICRU49NuclearStoppingModel : public G4VEmModel
+{
+
+public:
+
+  G4ICRU49NuclearStoppingModel(const G4String& nam = "ICRU49NucStopping");
+
+  virtual ~G4ICRU49NuclearStoppingModel();
+
+  virtual void Initialise(const G4ParticleDefinition*, 
+			  const G4DataVector&);
+
+  // main method to compute dEdx
+  virtual G4double ComputeDEDXPerVolume(const G4Material*,
+                                        const G4ParticleDefinition*,
+                                        G4double kineticEnergy,
+                                        G4double cutEnergy = DBL_MAX);
+
+  virtual void SampleSecondaries(std::vector<G4DynamicParticle*>*,
+				 const G4MaterialCutsCouple*,
+				 const G4DynamicParticle*, 
+				 G4double, G4double);
+
+  inline void SetFluctuationFlag(G4bool);
+
+private:
+
+  void InitialiseNuclearStopping();
+
+  G4double NuclearStoppingPower(G4double kineticEnergy,
+				G4double Z1, G4double Z2,
+				G4double A1, G4double A2);
+
+  //  hide assignment operator
+  G4ICRU49NuclearStoppingModel & operator=(const  G4ICRU49NuclearStoppingModel &right);
+  G4ICRU49NuclearStoppingModel(const  G4ICRU49NuclearStoppingModel&);
+
+  G4Pow* g4pow;
+
+  static G4double ed[104];
+  static G4double ad[104];
+
+  G4double theZieglerFactor;
+
+  // flags
+  G4bool   lossFlucFlag;
+};
+
+inline void G4ICRU49NuclearStoppingModel::SetFluctuationFlag(G4bool val)
+{
+  lossFlucFlag = val;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
+
Index: trunk/source/processes/electromagnetic/standard/include/G4LindhardPartition.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4LindhardPartition.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4LindhardPartition.hh	(revision 1197)
@@ -0,0 +1,69 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4LindhardPartition.hh,v 1.1 2009/11/05 11:42:24 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+
+/*
+ *  G4LindhardPartition.hh
+ *
+ *  Created by Marcus Mendenhall on 1/14/08.
+ *  Copyright 2008 Vanderbilt University, Nashville, TN, USA. All rights reserved.
+ *
+ */
+
+#include "globals.hh"
+
+class G4Material;
+
+class G4VNIELPartition 
+{
+public:
+	G4VNIELPartition() { }
+	virtual ~G4VNIELPartition() { }
+	
+	// return the fraction of the specified energy which will be deposited as NIEL
+	// if an incoming particle with z1, a1 is stopped in the specified material
+	// a1 is in atomic mass units, energy in native G4 energy units.
+	virtual G4double PartitionNIEL(
+		G4int z1, G4double a1, const G4Material *material, G4double energy
+	) const =0;
+};
+
+class G4LindhardRobinsonPartition : public G4VNIELPartition
+{
+public:
+	G4LindhardRobinsonPartition();
+	virtual ~G4LindhardRobinsonPartition() { }
+	
+	virtual G4double PartitionNIEL(
+		G4int z1, G4double a1, const G4Material *material, G4double energy
+	) const ;
+	
+	G4double z23[120];
+        size_t   max_z;
+};
+
Index: trunk/source/processes/electromagnetic/standard/include/G4NuclearStopping.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4NuclearStopping.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4NuclearStopping.hh	(revision 1197)
@@ -0,0 +1,109 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4NuclearStopping.hh,v 1.2 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -----------------------------------------------------------------------------
+//
+// GEANT4 Class header file
+//
+// Author:        Vladimir Ivanchenko
+//
+// Creation date: 20 July 2009
+// 
+// Modified:
+//
+//------------------------------------------------------------------------------
+//
+
+// class description
+//
+//  The class simulates nuclear stopping due to multiple scattering
+//
+// class description - end
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#ifndef G4NuclearStopping_h
+#define G4NuclearStopping_h 1
+
+#include "G4VEmProcess.hh"
+#include "globals.hh"
+#include "G4Track.hh"
+#include "G4Step.hh"
+#include "G4ParticleDefinition.hh"
+#include "G4GPILSelection.hh"
+#include "G4ParticleChangeForLoss.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4ICRU49NuclearStoppingModel;
+
+class G4NuclearStopping : public G4VEmProcess
+{
+
+public:    // with description
+
+  G4NuclearStopping(const G4String& processName="nuclearStopping");
+
+  virtual ~G4NuclearStopping();
+
+  // returns true for charged particles, false otherwise
+  G4bool IsApplicable (const G4ParticleDefinition& p);
+
+  // implementation of pure virtual method
+  G4double AlongStepGetPhysicalInteractionLength(const G4Track& track,
+						 G4double  previousStepSize,
+						 G4double  currentMinimumStep,
+						 G4double& proposedSafety,
+						 G4GPILSelection* selection);
+
+  // implementation of energy loss along step
+  G4VParticleChange* AlongStepDoIt(const G4Track& track,
+				   const G4Step&  step);
+
+  // Print few lines of informations about the process: validity range,
+  void PrintInfo();
+
+protected:
+
+  // This function initialise process
+  void InitialiseProcess(const G4ParticleDefinition*);
+
+private:       
+
+  G4ParticleChangeForLoss nParticleChange;
+
+  G4ICRU49NuclearStoppingModel* modelICRU49;
+
+  G4bool   isInitialized;
+
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#endif
Index: trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel92.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel92.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel92.hh	(revision 1197)
@@ -0,0 +1,230 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4UrbanMscModel92.hh,v 1.1 2009/11/01 13:04:12 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+//
+// GEANT4 Class header file
+//
+//
+// File name:     G4UrbanMscModel92
+//
+// Author:        Laszlo Urban
+//
+// Creation date: 06.03.2008
+//
+// Modifications:
+//
+// 28-10-2009  V.Ivanchenko moved G4UrbanMscModel to G4UrbanMscModel92, 
+//             now it is a frozen version of the Urban model corresponding 
+//             to g4 9.2
+//
+// Class Description:
+//
+// Implementation of the model of multiple scattering based on
+// H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
+
+// -------------------------------------------------------------------
+//
+
+#ifndef G4UrbanMscModel92_h
+#define G4UrbanMscModel92_h 1
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4VMscModel.hh"
+#include "G4PhysicsTable.hh"
+#include "G4MscStepLimitType.hh"
+
+class G4ParticleChangeForMSC;
+class G4SafetyHelper;
+class G4LossTableManager;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4UrbanMscModel92 : public G4VMscModel
+{
+
+public:
+
+  G4UrbanMscModel92(const G4String& nam = "UrbanMsc92");
+
+  virtual ~G4UrbanMscModel92();
+
+  void Initialise(const G4ParticleDefinition*, const G4DataVector&);
+
+  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
+				      G4double KineticEnergy,
+				      G4double AtomicNumber,
+				      G4double AtomicWeight=0., 
+				      G4double cut =0.,
+				      G4double emax=DBL_MAX);
+
+  void SampleScattering(const G4DynamicParticle*,
+			G4double safety);
+
+  G4double ComputeTruePathLengthLimit(const G4Track& track,
+				      G4PhysicsTable* theLambdaTable,
+				      G4double currentMinimalStep);
+
+  G4double ComputeGeomPathLength(G4double truePathLength);
+
+  G4double ComputeTrueStepLength(G4double geomStepLength);
+
+  G4double ComputeTheta0(G4double truePathLength,
+                         G4double KineticEnergy);
+
+private:
+
+  G4double SimpleScattering(G4double xmeanth, G4double x2meanth);
+
+  G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
+
+  G4double SampleDisplacement();
+
+  G4double LatCorrelation();
+
+  inline G4double GetLambda(G4double kinEnergy);
+
+  inline void SetParticle(const G4ParticleDefinition*);
+
+  inline void UpdateCache();
+
+  //  hide assignment operator
+  G4UrbanMscModel92 & operator=(const  G4UrbanMscModel92 &right);
+  G4UrbanMscModel92(const  G4UrbanMscModel92&);
+
+  const G4ParticleDefinition* particle;
+  G4ParticleChangeForMSC*     fParticleChange;
+
+  G4PhysicsTable*             theLambdaTable;
+  const G4MaterialCutsCouple* couple;
+  G4LossTableManager*         theManager;
+
+  G4double mass;
+  G4double charge,ChargeSquare;
+  G4double masslimite,lambdalimit,fr;
+
+  G4double taubig;
+  G4double tausmall;
+  G4double taulim;
+  G4double currentTau;
+  G4double tlimit;
+  G4double tlimitmin;
+  G4double tlimitminfix;
+  G4double tgeom;
+
+  G4double geombig;
+  G4double geommin;
+  G4double geomlimit;
+  G4double skindepth;
+  G4double smallstep;
+
+  G4double presafety;
+
+  G4double lambda0;
+  G4double lambdaeff;
+  G4double tPathLength;
+  G4double zPathLength;
+  G4double par1,par2,par3;
+
+  G4double stepmin;
+
+  G4double currentKinEnergy;
+  G4double currentRange; 
+  G4double rangeinit;
+  G4double currentRadLength;
+
+  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......
+
+inline
+G4double G4UrbanMscModel92::GetLambda(G4double e)
+{
+  G4double x;
+  if(theLambdaTable) {
+    G4bool b;
+    x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
+  } else {
+    x = CrossSection(couple,particle,e);
+  }
+  if(x > DBL_MIN) x = 1./x;
+  else            x = DBL_MAX;
+  return x;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline
+void G4UrbanMscModel92::SetParticle(const G4ParticleDefinition* p)
+{
+  if (p != particle) {
+    particle = p;
+    mass = p->GetPDGMass();
+    charge = p->GetPDGCharge()/eplus;
+    ChargeSquare = charge*charge;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline
+void G4UrbanMscModel92::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/G4UrbanMscModel93.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel93.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4UrbanMscModel93.hh	(revision 1197)
@@ -0,0 +1,235 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4UrbanMscModel93.hh,v 1.1 2009/11/01 13:04:12 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+//
+// GEANT4 Class header file
+//
+//
+// File name:     G4UrbanMscModel93
+//
+// Author:        Laszlo Urban
+//
+// Creation date: 06.03.2008
+//
+// Modifications:
+// 23-04-2009 L.Urban updated parameterization in UpdateCache method
+// 28-10-2009 V.Ivanchenko moved G4UrbanMscModel2 to G4UrbanMscModel93, 
+//            now it is a frozen version of the Urban model corresponding 
+//            to g4 9.3
+//
+//
+// Class Description:
+//
+// Implementation of the model of multiple scattering based on
+// H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
+
+// -------------------------------------------------------------------
+//
+
+#ifndef G4UrbanMscModel93_h
+#define G4UrbanMscModel93_h 1
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4VMscModel.hh"
+#include "G4PhysicsTable.hh"
+#include "G4MscStepLimitType.hh"
+
+class G4ParticleChangeForMSC;
+class G4SafetyHelper;
+class G4LossTableManager;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+class G4UrbanMscModel93 : public G4VMscModel
+{
+
+public:
+
+  G4UrbanMscModel93(const G4String& nam = "UrbanMsc93");
+
+  virtual ~G4UrbanMscModel93();
+
+  void Initialise(const G4ParticleDefinition*, const G4DataVector&);
+
+  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition* particle,
+				      G4double KineticEnergy,
+				      G4double AtomicNumber,
+				      G4double AtomicWeight=0., 
+				      G4double cut =0.,
+				      G4double emax=DBL_MAX);
+
+  void SampleScattering(const G4DynamicParticle*,
+			G4double safety);
+
+  G4double ComputeTruePathLengthLimit(const G4Track& track,
+				      G4PhysicsTable* theLambdaTable,
+				      G4double currentMinimalStep);
+
+  G4double ComputeGeomPathLength(G4double truePathLength);
+
+  G4double ComputeTrueStepLength(G4double geomStepLength);
+
+  G4double ComputeTheta0(G4double truePathLength,
+                         G4double KineticEnergy);
+
+private:
+
+  G4double SimpleScattering(G4double xmeanth, G4double x2meanth);
+
+  G4double SampleCosineTheta(G4double trueStepLength, G4double KineticEnergy);
+
+  G4double SampleDisplacement();
+
+  G4double LatCorrelation();
+
+  inline G4double GetLambda(G4double kinEnergy);
+
+  inline void SetParticle(const G4ParticleDefinition*);
+
+  inline void UpdateCache();
+
+  //  hide assignment operator
+  G4UrbanMscModel93 & operator=(const  G4UrbanMscModel93 &right);
+  G4UrbanMscModel93(const  G4UrbanMscModel93&);
+
+  const G4ParticleDefinition* particle;
+  G4ParticleChangeForMSC*     fParticleChange;
+
+  G4PhysicsTable*             theLambdaTable;
+  const G4MaterialCutsCouple* couple;
+  G4LossTableManager*         theManager;
+
+  G4double mass;
+  G4double charge,ChargeSquare;
+  G4double masslimite,lambdalimit,fr;
+
+  G4double taubig;
+  G4double tausmall;
+  G4double taulim;
+  G4double currentTau;
+  G4double tlimit;
+  G4double tlimitmin;
+  G4double tlimitminfix;
+  G4double tgeom;
+
+  G4double geombig;
+  G4double geommin;
+  G4double geomlimit;
+  G4double skindepth;
+  G4double smallstep;
+
+  G4double presafety;
+
+  G4double lambda0;
+  G4double lambdaeff;
+  G4double tPathLength;
+  G4double zPathLength;
+  G4double par1,par2,par3;
+
+  G4double stepmin;
+
+  G4double currentKinEnergy;
+  G4double currentRange; 
+  G4double rangeinit;
+  G4double currentRadLength;
+
+  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......
+
+inline
+G4double G4UrbanMscModel93::GetLambda(G4double e)
+{
+  G4double x;
+  if(theLambdaTable) {
+    G4bool b;
+    x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
+  } else {
+    x = CrossSection(couple,particle,e);
+  }
+  if(x > DBL_MIN) x = 1./x;
+  else            x = DBL_MAX;
+  return x;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline
+void G4UrbanMscModel93::SetParticle(const G4ParticleDefinition* p)
+{
+  if (p != particle) {
+    particle = p;
+    mass = p->GetPDGMass();
+    charge = p->GetPDGCharge()/eplus;
+    ChargeSquare = charge*charge;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+inline
+void G4UrbanMscModel93::UpdateCache()                                   
+{
+    lnZ = std::log(Zeff);
+    //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;
+                                              
+    Zold = Zeff;
+}
+
+#endif
+
Index: trunk/source/processes/electromagnetic/standard/include/G4alphaIonisation.hh
===================================================================
--- trunk/source/processes/electromagnetic/standard/include/G4alphaIonisation.hh	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/include/G4alphaIonisation.hh	(revision 1197)
@@ -0,0 +1,127 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4alphaIonisation.hh,v 1.1 2009/11/10 11:50:30 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class header file
+//
+//
+// File name:     G4alphaIonisation
+//
+// Author:        Vladimir Ivanchenko
+//
+// Creation date: 28.10.2009
+//
+// Modifications:
+//
+//
+// Class Description:
+//
+// This class manages the ionisation process for He ions. Effective charge,
+// nuclear stopping power, energy loss corrections are taken into account.
+// It inherites from G4VEnergyLossLoss.
+//
+
+// -------------------------------------------------------------------
+//
+
+#ifndef G4alphaIonisation_h
+#define G4alphaIonisation_h 1
+
+#include "G4VEnergyLossProcess.hh"
+
+class G4Material;
+
+class G4alphaIonisation : public G4VEnergyLossProcess
+{
+public:
+
+  G4alphaIonisation(const G4String& name = "alphaIoni");
+
+  virtual ~G4alphaIonisation();
+
+  virtual G4bool IsApplicable(const G4ParticleDefinition& p);
+
+  // Print out of the class parameters
+  virtual void PrintInfo();
+
+  void ActivateNuclearStopping(G4bool);
+
+protected:
+
+  virtual void InitialiseEnergyLossProcess(const G4ParticleDefinition*,
+					   const G4ParticleDefinition*);
+
+  virtual G4double MinPrimaryEnergy(const G4ParticleDefinition* p,
+				   const G4Material*, G4double cut);
+
+  inline G4double BetheBlochEnergyThreshold();
+
+  inline G4bool NuclearStoppingFlag();
+
+private:
+
+  // hide assignment operator
+  G4alphaIonisation & operator=(const G4alphaIonisation &right);
+  G4alphaIonisation(const G4alphaIonisation&);
+
+  const G4ParticleDefinition* theParticle;
+
+  G4double   mass;
+  G4double   ratio;
+  G4double   eth;
+
+  G4bool     isInitialised;
+  G4bool     nuclearStopping;
+};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline void G4alphaIonisation::ActivateNuclearStopping(G4bool val)
+{
+  nuclearStopping = val;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline G4double G4alphaIonisation::BetheBlochEnergyThreshold()
+{
+  return eth;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+inline G4bool G4alphaIonisation::NuclearStoppingFlag()
+{
+  return nuclearStopping;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#endif
Index: trunk/source/processes/electromagnetic/standard/src/G4ICRU49NuclearStoppingModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4ICRU49NuclearStoppingModel.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4ICRU49NuclearStoppingModel.cc	(revision 1197)
@@ -0,0 +1,311 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4ICRU49NuclearStoppingModel.cc,v 1.2 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:   G4ICRU49NuclearStoppingModel
+//
+// Author:      V.Ivanchenko 
+//
+// Creation date: 09.04.2008 from G4MuMscModel
+//
+// Modifications:
+//
+//
+// Class Description:
+//
+// Implementation of the model of ICRU'49 nuclear stopping
+
+// -------------------------------------------------------------------
+//
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4ICRU49NuclearStoppingModel.hh"
+#include "Randomize.hh"
+#include "G4LossTableManager.hh"
+#include "G4ParticleChangeForLoss.hh"
+#include "G4ElementVector.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4Step.hh"
+#include "G4Pow.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4ICRU49NuclearStoppingModel::ad[] = {0.0};
+G4double G4ICRU49NuclearStoppingModel::ed[] = {0.0};
+
+using namespace std;
+
+G4ICRU49NuclearStoppingModel::G4ICRU49NuclearStoppingModel(const G4String& nam) 
+  : G4VEmModel(nam),lossFlucFlag(false)
+{
+  theZieglerFactor = eV*cm2*1.0e-15;
+  g4pow = G4Pow::GetInstance();
+  if(ad[0] == 0.0) InitialiseNuclearStopping();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4ICRU49NuclearStoppingModel::~G4ICRU49NuclearStoppingModel()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4ICRU49NuclearStoppingModel::Initialise(const G4ParticleDefinition*, 
+					      const G4DataVector&)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void 
+G4ICRU49NuclearStoppingModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
+						const G4MaterialCutsCouple*,
+						const G4DynamicParticle*,
+						G4double, G4double)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double 
+G4ICRU49NuclearStoppingModel::ComputeDEDXPerVolume(const G4Material* mat,
+						   const G4ParticleDefinition* p,
+						   G4double kinEnergy,
+						   G4double)
+{
+  G4double nloss = 0.0;
+  if(kinEnergy <= 0.0) return nloss; 
+
+  // projectile
+  G4double m1 = p->GetPDGMass();
+  G4double z1 = std::fabs(p->GetPDGCharge()/eplus);
+
+  if(kinEnergy*proton_mass_c2/m1 > z1*z1*100.*MeV) return nloss;
+
+  // Projectile nucleus
+  m1 /= amu_c2;
+
+  //  loop for the elements in the material
+  G4int numberOfElements = mat->GetNumberOfElements();
+  const G4ElementVector* theElementVector = mat->GetElementVector();
+  const G4double* atomDensity  = mat->GetAtomicNumDensityVector();
+ 
+  for (G4int iel=0; iel<numberOfElements; iel++) {
+    const G4Element* element = (*theElementVector)[iel] ;
+    G4double z2 = element->GetZ();
+    G4double m2 = element->GetA()*mole/g ;
+    nloss += (NuclearStoppingPower(kinEnergy, z1, z2, m1, m2))
+           * atomDensity[iel] ;
+  }
+  nloss *= theZieglerFactor;
+  return nloss;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double 
+G4ICRU49NuclearStoppingModel::NuclearStoppingPower(G4double kineticEnergy,
+						   G4double z1, G4double z2,
+						   G4double m1, G4double m2)
+{
+  G4double energy = kineticEnergy/keV ;  // energy in keV
+  G4double nloss = 0.0;
+  G4double z12 = z1*z2;
+  G4int iz1 = G4int(z1);
+  G4int iz2 = G4int(z2);
+  
+  G4double rm;
+  if(iz1 > 1) rm = (m1 + m2) * ( g4pow->Z23(iz1) + g4pow->Z23(iz2) );
+  else        rm = (m1 + m2) * g4pow->Z13(iz2);
+
+  G4double er = 32.536 * m2 * energy / ( z12 * rm ) ;  // reduced energy
+
+  if (er >= ed[0]) nloss = ad[0];
+  else {
+    // the table is inverse in energy
+    for (G4int i=102; i>=0; --i)
+    {
+      if (er <= ed[i]) {
+	nloss = (ad[i] - ad[i+1])*(er - ed[i+1])/(ed[i] - ed[i+1]) + ad[i+1];
+	break;
+      }
+    }
+  }
+
+  // Stragling
+  if(lossFlucFlag) {
+    //    G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)*
+    //              (4.0 + 0.197*std::pow(er,-1.6991)+6.584*std::pow(er,-1.0494))) ;
+    G4double sig = 4.0 * m1 * m2 / ((m1 + m2)*(m1 + m2)*
+				    (4.0 + 0.197/(er*er) + 6.584/er));
+
+    nloss *= G4RandGauss::shoot(1.0,sig) ;
+    lossFlucFlag = false;
+  }
+   
+  nloss *= 8.462 * z12 * m1 / rm ; // Return to [ev/(10^15 atoms/cm^2]
+
+  if ( nloss < 0.0) nloss = 0.0 ;
+
+  return nloss;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4ICRU49NuclearStoppingModel::InitialiseNuclearStopping()
+{
+  const G4double nuca[104][2] = {
+  { 1.0E+8, 5.831E-8},
+  { 8.0E+7, 7.288E-8},
+  { 6.0E+7, 9.719E-8},
+  { 5.0E+7, 1.166E-7},
+  { 4.0E+7, 1.457E-7},
+  { 3.0E+7, 1.942E-7},
+  { 2.0E+7, 2.916E-7},
+  { 1.5E+7, 3.887E-7},
+
+  { 1.0E+7, 5.833E-7},
+  { 8.0E+6, 7.287E-7},
+  { 6.0E+6, 9.712E-7},
+  { 5.0E+6, 1.166E-6},
+  { 4.0E+6, 1.457E-6},
+  { 3.0E+6, 1.941E-6},
+  { 2.0E+6, 2.911E-6},
+  { 1.5E+6, 3.878E-6},
+
+  { 1.0E+6, 5.810E-6},
+  { 8.0E+5, 7.262E-6},
+  { 6.0E+5, 9.663E-6},
+  { 5.0E+5, 1.157E-5},
+  { 4.0E+5, 1.442E-5},
+  { 3.0E+5, 1.913E-5},
+  { 2.0E+5, 2.845E-5},
+  { 1.5E+5, 3.762E-5},
+
+  { 1.0E+5, 5.554E-5},
+  { 8.0E+4, 6.866E-5},
+  { 6.0E+4, 9.020E-5},
+  { 5.0E+4, 1.070E-4},
+  { 4.0E+4, 1.319E-4},
+  { 3.0E+4, 1.722E-4},
+  { 2.0E+4, 2.499E-4},
+  { 1.5E+4, 3.248E-4},
+
+  { 1.0E+4, 4.688E-4},
+  { 8.0E+3, 5.729E-4},
+  { 6.0E+3, 7.411E-4},
+  { 5.0E+3, 8.718E-4},
+  { 4.0E+3, 1.063E-3},
+  { 3.0E+3, 1.370E-3},
+  { 2.0E+3, 1.955E-3},
+  { 1.5E+3, 2.511E-3},
+
+  { 1.0E+3, 3.563E-3},
+  { 8.0E+2, 4.314E-3},
+  { 6.0E+2, 5.511E-3},
+  { 5.0E+2, 6.430E-3},
+  { 4.0E+2, 7.756E-3},
+  { 3.0E+2, 9.855E-3},
+  { 2.0E+2, 1.375E-2},
+  { 1.5E+2, 1.736E-2},
+
+  { 1.0E+2, 2.395E-2},
+  { 8.0E+1, 2.850E-2},
+  { 6.0E+1, 3.552E-2},
+  { 5.0E+1, 4.073E-2},
+  { 4.0E+1, 4.802E-2},
+  { 3.0E+1, 5.904E-2},
+  { 1.5E+1, 9.426E-2},
+
+  { 1.0E+1, 1.210E-1},
+  { 8.0E+0, 1.377E-1},
+  { 6.0E+0, 1.611E-1},
+  { 5.0E+0, 1.768E-1},
+  { 4.0E+0, 1.968E-1},
+  { 3.0E+0, 2.235E-1},
+  { 2.0E+0, 2.613E-1},
+  { 1.5E+0, 2.871E-1},
+
+  { 1.0E+0, 3.199E-1},
+  { 8.0E-1, 3.354E-1},
+  { 6.0E-1, 3.523E-1},
+  { 5.0E-1, 3.609E-1},
+  { 4.0E-1, 3.693E-1},
+  { 3.0E-1, 3.766E-1},
+  { 2.0E-1, 3.803E-1},
+  { 1.5E-1, 3.788E-1},
+
+  { 1.0E-1, 3.711E-1},
+  { 8.0E-2, 3.644E-1},
+  { 6.0E-2, 3.530E-1},
+  { 5.0E-2, 3.444E-1},
+  { 4.0E-2, 3.323E-1},
+  { 3.0E-2, 3.144E-1},
+  { 2.0E-2, 2.854E-1},
+  { 1.5E-2, 2.629E-1},
+
+  { 1.0E-2, 2.298E-1},
+  { 8.0E-3, 2.115E-1},
+  { 6.0E-3, 1.883E-1},
+  { 5.0E-3, 1.741E-1},
+  { 4.0E-3, 1.574E-1},
+  { 3.0E-3, 1.372E-1},
+  { 2.0E-3, 1.116E-1},
+  { 1.5E-3, 9.559E-2},
+
+  { 1.0E-3, 7.601E-2},
+  { 8.0E-4, 6.668E-2},
+  { 6.0E-4, 5.605E-2},
+  { 5.0E-4, 5.008E-2},
+  { 4.0E-4, 4.352E-2},
+  { 3.0E-4, 3.617E-2},
+  { 2.0E-4, 2.768E-2},
+  { 1.5E-4, 2.279E-2},
+
+  { 1.0E-4, 1.723E-2},
+  { 8.0E-5, 1.473E-2},
+  { 6.0E-5, 1.200E-2},
+  { 5.0E-5, 1.052E-2},
+  { 4.0E-5, 8.950E-3},
+  { 3.0E-5, 7.246E-3},
+  { 2.0E-5, 5.358E-3},
+  { 1.5E-5, 4.313E-3},
+  { 0.0, 3.166E-3}
+  };
+
+  for(G4int i=0; i<104; ++i) {
+    ed[i] = nuca[i][0];
+    ad[i] = nuca[i][1];
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4LindhardPartition.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4LindhardPartition.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4LindhardPartition.cc	(revision 1197)
@@ -0,0 +1,100 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4LindhardPartition.cc,v 1.1 2009/11/05 11:42:24 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+
+/*
+ *  G4LindhardPartition.cc
+ *
+ *  Created by Marcus Mendenhall on 1/14/08.
+ *  Copyright 2008 Vanderbilt University, Nashville, TN, USA. All rights reserved.
+ *
+ */
+
+#include "G4LindhardPartition.hh"
+#include "G4Material.hh"
+#include "G4Element.hh"
+
+/*
+for a first cut, we will compute NIEL from a Lindhard-Robinson partition
+based on the most abundant element in the material.
+
+this is from IEEE Trans. Nucl Science Vol. 48 No.1 February 2001 page 162++
+Insoo Jun, "Effects of Secondary Particles on the Total Dose..."
+and, by reference,
+Lindhard, Nielsen, Scharff & Thompson, 
+"Integral Equations Governing Radiation Efects...", 
+Mat. Fys. Medd. Dan. Vid. Selsk. vol 33 #10, pp1-42, 1963
+and
+Robinson, "The dependence of radiation effects on primary recoil energy",
+in Proc. Int. Conf. Radiation-Induced Voids in Metal,  
+Albany, NY 1972 pp. 397-439
+def lindhard_robinson(z1, a1, z2, a2, ke):
+el=30.724*z1*z2*math.sqrt(z1**0.6667+z2**0.6667)*(a1+a2)/a2
+fl=0.0793*z1**0.6667*math.sqrt(z2)*(a1+a2)**1.5/
+((z1**0.6667+z2**0.6667)**0.75*a1**1.5*math.sqrt(a2))
+eps=ke*(1.0/el)
+return 1.0/(1+fl*(3.4008*eps**0.16667+0.40244*eps**0.75+eps))
+*/
+
+G4LindhardRobinsonPartition::G4LindhardRobinsonPartition() 
+{
+  max_z = 120;
+  for(size_t i=1; i<max_z; i++) {z23[i]=std::pow((G4double)i, 2./3.);}
+}
+
+G4double G4LindhardRobinsonPartition::PartitionNIEL(
+   G4int z1, G4double a1, const G4Material *material, G4double energy) const
+{
+  size_t nMatElements = material->GetNumberOfElements();
+	
+  const G4double *atomDensities=material->GetVecNbOfAtomsPerVolume();
+  G4double maxdens=0.0;
+  size_t maxindex=0;
+  for (size_t k=0 ; k < nMatElements ; k++ )
+    {
+      if(atomDensities[k] > maxdens) {
+	maxdens=atomDensities[k];
+	maxindex=k;
+      }
+    }
+  const G4Element *element=material->GetElement(maxindex);
+
+  G4int z2=G4int(element->GetZ());
+	
+  G4double a2=element->GetA()/(Avogadro*amu);
+	
+  G4double zpow=z23[z1]+z23[z2];
+  G4double asum=a1+a2;
+	
+  G4double el=30.724*z1*z2*std::sqrt(zpow)*asum/a2;
+  G4double fl=0.0793*z23[z1]*std::sqrt(z2*asum*asum*asum/(a1*a1*a1*a2))/std::pow(zpow, 0.75);
+  G4double eps=(energy/eV)*(1.0/el);
+
+  return 1.0/(1+fl*(3.4008*std::pow(eps, 0.16667)+0.40244*std::pow(eps, 0.75)+eps));
+}
+
Index: trunk/source/processes/electromagnetic/standard/src/G4NuclearStopping.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4NuclearStopping.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4NuclearStopping.cc	(revision 1197)
@@ -0,0 +1,142 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4NuclearStopping.cc,v 1.2 2009/11/10 19:25:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -----------------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+// File name:     G4NuclearStopping
+//
+// Author:        Vladimir Ivanchenko
+//
+// Creation date: 20 July 2009
+// 
+// Modified:
+//
+// -----------------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4NuclearStopping.hh"
+#include "G4ICRU49NuclearStoppingModel.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+using namespace std;
+
+G4NuclearStopping::G4NuclearStopping(const G4String& processName)
+  : G4VEmProcess(processName)
+{
+  isInitialized = false;  
+  SetProcessSubType(fCoulombScattering);
+  enableAlongStepDoIt = true; 
+  modelICRU49 = 0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4NuclearStopping::~G4NuclearStopping()
+{
+  delete modelICRU49;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4bool G4NuclearStopping::IsApplicable (const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && !p.IsShortLived());
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4NuclearStopping::InitialiseProcess(const G4ParticleDefinition*)
+{
+  if(!isInitialized) {
+    isInitialized = true;
+    SetBuildTableFlag(false);
+    //    SetSecondaryParticle(G4Electron::Electron());
+    if(!Model()) {
+      modelICRU49 = new G4ICRU49NuclearStoppingModel();
+      SetModel(modelICRU49);
+    }
+    AddEmModel(1, Model());
+    Model()->SetParticleChange(&nParticleChange);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double 
+G4NuclearStopping::AlongStepGetPhysicalInteractionLength(
+     const G4Track&, G4double, G4double, G4double&, G4GPILSelection* selection)
+{
+  *selection = NotCandidateForSelection;
+  return DBL_MAX;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4VParticleChange* G4NuclearStopping::AlongStepDoIt(const G4Track& track,
+						    const G4Step&  step)
+{
+  nParticleChange.InitializeForAlongStep(track);
+  G4double length = step.GetStepLength();
+  if(length > 0.0) {
+
+    G4double T1 = step.GetPreStepPoint()->GetKineticEnergy();
+    G4double T2 = step.GetPostStepPoint()->GetKineticEnergy();
+    if(T2 > 0.0) {
+
+      // primary
+      G4double T = 0.5*(T1 + T2);
+      const G4ParticleDefinition* part = track.GetDefinition();
+      const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 
+      G4VEmModel* mod = SelectModel(T, couple->GetIndex());
+
+      // sample stopping
+      if(modelICRU49) modelICRU49->SetFluctuationFlag(true);
+      G4double nloss = 
+	length*mod->ComputeDEDXPerVolume(couple->GetMaterial(), part, T);
+      if(nloss > T1) nloss = T1;
+      nParticleChange.SetProposedKineticEnergy(T1 - nloss);
+      nParticleChange.ProposeLocalEnergyDeposit(nloss);
+      nParticleChange.ProposeNonIonizingEnergyDeposit(nloss);
+      if(modelICRU49) modelICRU49->SetFluctuationFlag(false);
+    }
+  }
+  return &nParticleChange;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4NuclearStopping::PrintInfo()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel92.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel92.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel92.cc	(revision 1197)
@@ -0,0 +1,1056 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//
+// $Id: G4UrbanMscModel92.cc,v 1.1 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:   G4UrbanMscModel92
+//
+// Author:      Laszlo Urban
+//
+// Creation date: 06.03.2008
+//
+// Modifications:
+//
+// 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-10-09  V.Ivanchenko moved G4UrbanMscModel to G4UrbanMscModel92, 
+//           now it is a frozen version of the Urban model corresponding 
+//           to g4 9.2
+
+// Class Description:
+//
+// Implementation of the model of multiple scattering based on
+// H.W.Lewis Phys Rev 78 (1950) 526 and others
+
+// -------------------------------------------------------------------
+// In its present form the model can be  used for simulation 
+//   of the e-/e+, muon and charged hadron multiple scattering
+//
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4UrbanMscModel92.hh"
+#include "Randomize.hh"
+#include "G4Electron.hh"
+#include "G4LossTableManager.hh"
+#include "G4ParticleChangeForMSC.hh"
+
+#include "G4Poisson.hh"
+#include "globals.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+using namespace std;
+
+G4UrbanMscModel92::G4UrbanMscModel92(const G4String& nam)
+  : G4VMscModel(nam),
+    isInitialized(false)
+{
+  masslimite    = 0.6*MeV;
+  lambdalimit   = 1.*mm;
+  fr            = 0.02;
+  //  facsafety     = 0.3;
+  taubig        = 8.0;
+  tausmall      = 1.e-16;
+  taulim        = 1.e-6;
+  currentTau    = taulim;
+  tlimitminfix  = 1.e-6*mm;            
+  stepmin       = tlimitminfix;
+  smallstep     = 1.e10;
+  currentRange  = 0. ;
+  rangeinit     = 0.;
+  tlimit        = 1.e10*mm;
+  tlimitmin     = 10.*tlimitminfix;            
+  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(); 
+  inside        = false;  
+  insideskin    = false;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4UrbanMscModel92::~G4UrbanMscModel92()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4UrbanMscModel92::Initialise(const G4ParticleDefinition* p,
+				   const G4DataVector&)
+{
+  skindepth = skin*stepmin;
+  if(isInitialized) return;
+  // set values of some data members
+  SetParticle(p);
+
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
+
+  isInitialized = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::ComputeCrossSectionPerAtom( 
+                             const G4ParticleDefinition* part,
+                                   G4double KineticEnergy,
+                                   G4double AtomicNumber,G4double,
+				   G4double, G4double)
+{
+  const G4double sigmafactor = twopi*classic_electr_radius*classic_electr_radius;
+  const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
+                            Bohr_radius*Bohr_radius/(hbarc*hbarc);
+  const G4double epsmin = 1.e-4 , epsmax = 1.e10;
+
+  const G4double Zdat[15] = { 4.,  6., 13., 20., 26., 29., 32., 38., 47.,
+                             50., 56., 64., 74., 79., 82. };
+
+  const G4double Tdat[22] = { 100*eV,  200*eV,  400*eV,  700*eV,
+                               1*keV,   2*keV,   4*keV,   7*keV,
+			      10*keV,  20*keV,  40*keV,  70*keV,
+                             100*keV, 200*keV, 400*keV, 700*keV,
+                               1*MeV,   2*MeV,   4*MeV,   7*MeV,
+			      10*MeV,  20*MeV};
+
+  // corr. factors for e-/e+ lambda for T <= Tlim
+          G4double celectron[15][22] =
+          {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
+            1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
+            1.112,1.108,1.100,1.093,1.089,1.087            },
+           {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
+            1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
+            1.109,1.105,1.097,1.090,1.086,1.082            },
+           {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
+            1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
+            1.131,1.124,1.113,1.104,1.099,1.098            },
+           {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
+            1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
+            1.112,1.105,1.096,1.089,1.085,1.098            },
+           {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
+            1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
+            1.073,1.070,1.064,1.059,1.056,1.056            },
+           {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
+            1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
+            1.074,1.070,1.063,1.059,1.056,1.052            },
+           {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
+            1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
+            1.068,1.064,1.059,1.054,1.051,1.050            },
+           {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
+            1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
+            1.039,1.037,1.034,1.031,1.030,1.036            },
+           {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
+            1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
+            1.031,1.028,1.024,1.022,1.021,1.024            },
+           {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
+            1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
+            1.020,1.017,1.015,1.013,1.013,1.020            },
+           {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
+            1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
+            0.995,0.993,0.993,0.993,0.993,1.011            },
+           {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
+            1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
+            0.974,0.972,0.973,0.974,0.975,0.987            },
+           {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
+            1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
+            0.950,0.947,0.949,0.952,0.954,0.963            },
+           {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
+            1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
+            0.941,0.938,0.940,0.944,0.946,0.954            },
+           {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
+            1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
+            0.933,0.930,0.933,0.936,0.939,0.949            }};
+	    
+           G4double cpositron[15][22] = {
+           {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
+            1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
+            1.131,1.126,1.117,1.108,1.103,1.100            },
+           {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
+            1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
+            1.138,1.132,1.122,1.113,1.108,1.102            },
+           {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
+            1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
+            1.203,1.190,1.173,1.159,1.151,1.145            },
+           {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
+            1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
+            1.225,1.210,1.191,1.175,1.166,1.174            },
+           {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
+            1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
+            1.217,1.203,1.184,1.169,1.160,1.151            },
+           {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
+            1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
+            1.237,1.222,1.201,1.184,1.174,1.159            },
+           {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
+            1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
+            1.252,1.234,1.212,1.194,1.183,1.170            },
+           {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
+            2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
+            1.254,1.237,1.214,1.195,1.185,1.179            },
+           {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
+            2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
+            1.312,1.288,1.258,1.235,1.221,1.205            },
+           {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
+            2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
+            1.320,1.294,1.264,1.240,1.226,1.214            },
+           {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
+            2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
+            1.328,1.302,1.270,1.245,1.231,1.233            },
+           {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
+            2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
+            1.361,1.330,1.294,1.267,1.251,1.239            },
+           {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
+            3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
+            1.409,1.372,1.330,1.298,1.280,1.258            },
+           {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
+            3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
+            1.442,1.400,1.354,1.319,1.299,1.272            },
+           {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
+            3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
+            1.456,1.412,1.364,1.328,1.307,1.282            }};
+
+  //data/corrections for T > Tlim  
+  G4double Tlim = 10.*MeV;
+  G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
+                      ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
+  G4double bg2lim   = Tlim*(Tlim+2.*electron_mass_c2)/
+                      (electron_mass_c2*electron_mass_c2);
+
+  G4double sig0[15] = {0.2672*barn,  0.5922*barn, 2.653*barn,  6.235*barn,
+                      11.69*barn  , 13.24*barn  , 16.12*barn, 23.00*barn ,
+		      35.13*barn  , 39.95*barn  , 50.85*barn, 67.19*barn ,
+                      91.15*barn  , 104.4*barn  , 113.1*barn};
+		      		       
+  G4double hecorr[15] = {120.70, 117.50, 105.00, 92.92, 79.23,  74.510,  68.29,
+                          57.39,  41.97,  36.14, 24.53, 10.21,  -7.855, -16.84,
+			 -22.30};
+
+  G4double sigma;
+  SetParticle(part);
+
+  G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23);
+
+  // correction if particle .ne. e-/e+
+  // compute equivalent kinetic energy
+  // lambda depends on p*beta ....
+
+  G4double eKineticEnergy = KineticEnergy;
+
+  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 eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
+  G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
+                                 /(eTotalEnergy*eTotalEnergy);
+  G4double bg2   = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
+                                 /(electron_mass_c2*electron_mass_c2);
+
+  G4double eps = epsfactor*bg2/Z23;
+
+  if     (eps<epsmin)  sigma = 2.*eps*eps;
+  else if(eps<epsmax)  sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
+  else                 sigma = log(2.*eps)-1.+1./eps;
+
+  sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
+
+  // interpolate in AtomicNumber and beta2 
+  G4double c1,c2,cc1,cc2,corr;
+
+  // get bin number in Z
+  G4int iZ = 14;
+  while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
+  if (iZ==14)                               iZ = 13;
+  if (iZ==-1)                               iZ = 0 ;
+
+  G4double Z1 = Zdat[iZ];
+  G4double Z2 = Zdat[iZ+1];
+  G4double ratZ = (AtomicNumber-Z1)*(AtomicNumber+Z1)/
+                  ((Z2-Z1)*(Z2+Z1));
+
+  if(eKineticEnergy <= Tlim) 
+  {
+    // get bin number in T (beta2)
+    G4int iT = 21;
+    while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
+    if(iT==21)                                  iT = 20;
+    if(iT==-1)                                  iT = 0 ;
+
+    //  calculate betasquare values
+    G4double T = Tdat[iT],   E = T + electron_mass_c2;
+    G4double b2small = T*(E+electron_mass_c2)/(E*E);
+
+    T = Tdat[iT+1]; E = T + electron_mass_c2;
+    G4double b2big = T*(E+electron_mass_c2)/(E*E);
+    G4double ratb2 = (beta2-b2small)/(b2big-b2small);
+
+    if (charge < 0.)
+    {
+       c1 = celectron[iZ][iT];
+       c2 = celectron[iZ+1][iT];
+       cc1 = c1+ratZ*(c2-c1);
+
+       c1 = celectron[iZ][iT+1];
+       c2 = celectron[iZ+1][iT+1];
+       cc2 = c1+ratZ*(c2-c1);
+
+       corr = cc1+ratb2*(cc2-cc1);
+
+       sigma *= sigmafactor/corr;
+    }
+    else              
+    {
+       c1 = cpositron[iZ][iT];
+       c2 = cpositron[iZ+1][iT];
+       cc1 = c1+ratZ*(c2-c1);
+
+       c1 = cpositron[iZ][iT+1];
+       c2 = cpositron[iZ+1][iT+1];
+       cc2 = c1+ratZ*(c2-c1);
+
+       corr = cc1+ratb2*(cc2-cc1);
+
+       sigma *= sigmafactor/corr;
+    }
+  }
+  else
+  {
+    c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
+    c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
+    if((AtomicNumber >= Z1) && (AtomicNumber <= Z2))
+      sigma = c1+ratZ*(c2-c1) ;
+    else if(AtomicNumber < Z1)
+      sigma = AtomicNumber*AtomicNumber*c1/(Z1*Z1);
+    else if(AtomicNumber > Z2)
+      sigma = AtomicNumber*AtomicNumber*c2/(Z2*Z2);
+  }
+  return sigma;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::ComputeTruePathLengthLimit(
+                             const G4Track& track,
+			     G4PhysicsTable* theTable,
+			     G4double currentMinimalStep)
+{
+  tPathLength = currentMinimalStep;
+  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
+  G4StepStatus stepStatus = sp->GetStepStatus();
+
+  const G4DynamicParticle* dp = track.GetDynamicParticle();
+
+  if(stepStatus == fUndefined) {
+    inside = false;
+    insideskin = false;
+    tlimit = geombig;
+    SetParticle( dp->GetDefinition() );
+  }
+
+  theLambdaTable = theTable;
+  couple = track.GetMaterialCutsCouple();
+  currentMaterialIndex = couple->GetIndex();
+  currentKinEnergy = dp->GetKineticEnergy();
+  currentRange = 
+    theManager->GetRangeFromRestricteDEDX(particle,currentKinEnergy,couple);
+  lambda0 = GetLambda(currentKinEnergy);
+
+  // stop here if small range particle
+  if(inside) return tPathLength;            
+  
+  if(tPathLength > currentRange) tPathLength = currentRange;
+
+  presafety = sp->GetSafety();
+
+  //    G4cout << "G4UrbanMscModel92::ComputeTruePathLengthLimit tPathLength= " 
+  //	   <<tPathLength<<" safety= " << presafety
+  //	   << " range= " <<currentRange<<G4endl;
+
+  // far from geometry boundary
+  if(currentRange < presafety)
+    {
+      inside = true;
+      return tPathLength;  
+    }
+
+  // standard  version
+  //
+  if (steppingAlgorithm == fUseDistanceToBoundary)
+    {
+      //compute geomlimit and presafety 
+      G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
+
+      // is it far from boundary ?
+      if(currentRange < presafety)
+	{
+	  inside = true;
+	  return tPathLength;   
+	}
+
+      smallstep += 1.;
+      insideskin = false;
+
+      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;
+
+          //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  
+      //     << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
+
+      // shortcut
+      if((tPathLength < tlimit) && (tPathLength < presafety) &&
+         (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
+	return tPathLength;   
+
+      // step reduction near to boundary
+      if(smallstep < skin)
+	{
+	  tlimit = stepmin;
+	  insideskin = true;
+	}
+      else if(geomlimit < geombig)
+	{
+	  if(geomlimit > skindepth)
+	    {
+	      if(tlimit > geomlimit-0.999*skindepth)
+		tlimit = geomlimit-0.999*skindepth;
+	    }
+	  else
+	    {
+	      insideskin = true;
+	      if(tlimit > stepmin) tlimit = stepmin;
+	    }
+	}
+
+      if(tlimit < stepmin) tlimit = stepmin;
+
+      if(tPathLength > tlimit) tPathLength = tlimit  ; 
+
+    }
+    // for 'normal' simulation with or without magnetic field 
+    //  there no small step/single scattering at boundaries
+  else if(steppingAlgorithm == fUseSafety)
+    {
+      // compute presafety again if presafety <= 0 and no boundary
+      // i.e. when it is needed for optimization purposes
+      if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
+	presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
+
+      // is far from boundary
+      if(currentRange < presafety)
+        {
+          inside = true;
+          return tPathLength;  
+        }
+
+      if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
+      {
+        rangeinit = currentRange;
+        fr = facrange;
+        // 9.1 like stepping for e+/e- only (not for muons,hadrons)
+        if(mass < masslimite) 
+        {
+          if(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;
+    }
+  
+  // version similar to 7.1 (needed for some experiments)
+  else
+    {
+      if (stepStatus == fGeomBoundary)
+	{
+	  if (currentRange > lambda0) tlimit = facrange*currentRange;
+	  else                        tlimit = facrange*lambda0;
+
+	  if(tlimit < tlimitmin) tlimit = tlimitmin;
+	  if(tPathLength > tlimit) tPathLength = tlimit;
+	}
+    }
+  // G4cout << "tPathLength= " << tPathLength << "  geomlimit= " << geomlimit 
+  //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
+  return tPathLength ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::ComputeGeomPathLength(G4double)
+{
+  lambdaeff = lambda0;
+  par1 = -1. ;  
+  par2 = par3 = 0. ;  
+
+  //  do the true -> geom transformation
+  zPathLength = tPathLength;
+
+  // z = t for very small tPathLength
+  if(tPathLength < tlimitminfix) return zPathLength;
+
+  // this correction needed to run MSC with eIoni and eBrem inactivated
+  // and makes no harm for a normal run
+  if(tPathLength > currentRange)
+    tPathLength = currentRange ;
+
+  G4double tau   = tPathLength/lambda0 ;
+
+  if ((tau <= tausmall) || insideskin) {
+    zPathLength  = tPathLength;
+    if(zPathLength > lambda0) zPathLength = lambda0;
+    return zPathLength;
+  }
+
+  G4double zmean = tPathLength;
+  if (tPathLength < currentRange*dtrl) {
+    if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
+    else             zmean = lambda0*(1.-exp(-tau));
+  } else if(currentKinEnergy < mass)  {
+    par1 = 1./currentRange ;
+    par2 = 1./(par1*lambda0) ;
+    par3 = 1.+par2 ;
+    if(tPathLength < currentRange)
+      zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
+    else
+      zmean = 1./(par1*par3) ;
+  } else {
+    G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,couple);
+    G4double lambda1 = GetLambda(T1);
+
+    par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
+    par2 = 1./(par1*lambda0) ;
+    par3 = 1.+par2 ;
+    zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
+  }
+
+  zPathLength = zmean ;
+
+  //  sample z
+  if(samplez)
+  {
+    const G4double  ztmax = 0.99 ;
+    G4double zt = zmean/tPathLength ;
+
+    if (tPathLength > stepmin && zt < ztmax)              
+    {
+      G4double u,cz1;
+      if(zt >= third)
+      {
+        G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
+        cz1 = 1.+cz ;
+        G4double u0 = cz/cz1 ;
+        G4double grej ;
+        do {
+            u = exp(log(G4UniformRand())/cz1) ;
+            grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
+           } while (grej < G4UniformRand()) ;
+      }
+      else
+      {
+        cz1 = 1./zt-1.;
+        u = 1.-exp(log(G4UniformRand())/cz1) ;
+      }
+      zPathLength = tPathLength*u ;
+    }
+  }
+
+  if(zPathLength > lambda0) zPathLength = lambda0;
+
+  return zPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::ComputeTrueStepLength(G4double geomStepLength)
+{
+  // step defined other than transportation 
+  if(geomStepLength == zPathLength && tPathLength <= currentRange)
+    return tPathLength;
+
+  // t = z for very small step
+  zPathLength = geomStepLength;
+  tPathLength = geomStepLength;
+  if(geomStepLength < tlimitminfix) return tPathLength;
+  
+  // recalculation
+  if((geomStepLength > lambda0*tausmall) && !insideskin)
+  {
+    if(par1 <  0.)
+      tPathLength = -lambda0*log(1.-geomStepLength/lambda0) ;
+    else 
+    {
+      if(par1*par3*geomStepLength < 1.)
+        tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
+      else 
+        tPathLength = currentRange;
+    }  
+  }
+  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
+
+  return tPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::ComputeTheta0(G4double trueStepLength,
+					  G4double KineticEnergy)
+{
+  // for all particles take the width of the central part
+  //  from a  parametrization similar to the Highland formula
+  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
+  const G4double c_highland = 13.6*MeV ;
+  G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
+                         KineticEnergy*(KineticEnergy+2.*mass)/
+                      ((currentKinEnergy+mass)*(KineticEnergy+mass)));
+  y = trueStepLength/currentRadLength;
+  G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
+  y = log(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;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4UrbanMscModel92::SampleScattering(const G4DynamicParticle* dynParticle,
+					 G4double safety)
+{
+  G4double kineticEnergy = dynParticle->GetKineticEnergy();
+
+  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) ||
+     (tPathLength/tausmall < lambda0)) return;
+
+  G4double cth  = SampleCosineTheta(tPathLength,kineticEnergy);
+
+  // protection against 'bad' cth values
+  if(std::abs(cth) > 1.) return;
+
+  G4double sth  = sqrt((1.0 - cth)*(1.0 + cth));
+  G4double phi  = twopi*G4UniformRand();
+  G4double dirx = sth*cos(phi);
+  G4double diry = sth*sin(phi);
+
+  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
+  G4ThreeVector newDirection(dirx,diry,cth);
+  newDirection.rotateUz(oldDirection);
+  fParticleChange->ProposeMomentumDirection(newDirection);
+
+  if (latDisplasment && safety > tlimitminfix) {
+
+    G4double r = SampleDisplacement();
+    /*
+    G4cout << "G4UrbanMscModel92::SampleSecondaries: e(MeV)= " << kineticEnergy
+	   << " sinTheta= " << sth << " r(mm)= " << r
+           << " trueStep(mm)= " << tPathLength
+           << " geomStep(mm)= " << zPathLength
+           << G4endl;
+    */
+    if(r > 0.)
+      {
+        G4double latcorr = LatCorrelation();
+        if(latcorr > r) latcorr = r;
+
+        // sample direction of lateral displacement
+        // compute it from the lateral correlation
+        G4double Phi = 0.;
+        if(std::abs(r*sth) < latcorr)
+          Phi  = twopi*G4UniformRand();
+        else
+        {
+          G4double psi = std::acos(latcorr/(r*sth));
+          if(G4UniformRand() < 0.5)
+            Phi = phi+psi;
+          else
+            Phi = phi-psi;
+        }
+
+        dirx = std::cos(Phi);
+        diry = std::sin(Phi);
+
+        G4ThreeVector latDirection(dirx,diry,0.0);
+        latDirection.rotateUz(oldDirection);
+
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::SampleCosineTheta(G4double trueStepLength,
+					      G4double KineticEnergy)
+{
+  G4double cth = 1. ;
+  G4double tau = trueStepLength/lambda0 ;
+
+  Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
+         couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
+
+  if(Zold != Zeff)  
+    UpdateCache();
+
+  if(insideskin)
+  {
+    //no scattering, single or plural scattering
+    G4double mean = trueStepLength/stepmin ;
+
+    G4int n = G4Poisson(mean);
+    if(n > 0)
+    {
+      //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;
+      G4double sx=0.,sy=0.,sz=0.;
+      for(G4int i=1; i<=n; i++)
+      {
+        ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
+        if(ct < -1.) ct = -1.;
+        if(ct >  1.) ct =  1.; 
+        st = sqrt(1.-ct*ct);
+        phi = twopi*G4UniformRand();
+        sx += st*cos(phi);
+        sy += st*sin(phi);
+        sz += ct;
+      }
+      cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
+    }
+  }
+  else
+  {
+    if(trueStepLength >= currentRange*dtrl) 
+    {
+      if(par1*trueStepLength < 1.)
+	tau = -par2*log(1.-par1*trueStepLength) ;
+      // for the case if ioni/brems are inactivated
+      // see the corresponding condition in ComputeGeomPathLength 
+      else if(1.-KineticEnergy/currentKinEnergy > taulim)
+	tau = taubig ;
+    }
+    currentTau = tau ;
+    lambdaeff = trueStepLength/currentTau;
+    currentRadLength = couple->GetMaterial()->GetRadlen();
+
+    if (tau >= taubig) cth = -1.+2.*G4UniformRand();
+    else if (tau >= tausmall)
+    {
+      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);
+
+      // 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(-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)
+      {
+        if(G4UniformRand() < prob)
+          cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
+        else
+          cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
+      }
+      else
+        cth = -1.+2.*G4UniformRand();
+    }
+  }  
+  return cth ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::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;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::SampleDisplacement()
+{
+  const G4double kappa = 2.5;
+  const G4double kappapl1 = kappa+1.;
+  const G4double kappami1 = kappa-1.;
+  G4double rmean = 0.0;
+  if ((currentTau >= tausmall) && !insideskin) {
+    if (currentTau < taulim) {
+      rmean = kappa*currentTau*currentTau*currentTau*
+             (1.-kappapl1*currentTau*0.25)/6. ;
+
+    } else {
+      G4double etau = 0.0;
+      if (currentTau<taubig) etau = exp(-currentTau);
+      rmean = -kappa*currentTau;
+      rmean = -exp(rmean)/(kappa*kappami1);
+      rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
+    }
+    if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean/3.0);
+    else          rmean = 0.;
+  }
+
+  // protection against z > t ...........................
+  if(rmean > 0.) {
+    G4double zt = (tPathLength-zPathLength)*(tPathLength+zPathLength);
+    if(zt <= 0.)
+      rmean = 0.;
+    else if(rmean*rmean > zt)
+      rmean = sqrt(zt);
+  }
+  return rmean;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel92::LatCorrelation()
+{
+  const G4double kappa = 2.5;
+  const G4double kappami1 = kappa-1.;
+
+  G4double latcorr = 0.;
+  if((currentTau >= tausmall) && !insideskin)
+  {
+    if(currentTau < taulim)
+      latcorr = lambdaeff*kappa*currentTau*currentTau*
+                (1.-(kappa+1.)*currentTau/3.)/3.;
+    else
+    {
+      G4double etau = 0.;
+      if(currentTau < taubig) etau = exp(-currentTau);
+      latcorr = -kappa*currentTau;
+      latcorr = exp(latcorr)/kappami1;
+      latcorr += 1.-kappa*etau/kappami1 ;
+      latcorr *= 2.*lambdaeff/3. ;
+    }
+  }
+
+  return latcorr;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel93.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel93.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4UrbanMscModel93.cc	(revision 1197)
@@ -0,0 +1,1062 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//
+// $Id: G4UrbanMscModel93.cc,v 1.1 2009/11/01 13:05:01 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:   G4UrbanMscModel93
+//
+// Author:      Laszlo Urban
+//
+// Creation date: 06.03.2008
+//
+// Modifications:
+//
+// 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-10-09  V.Ivanchenko moved G4UrbanMscModel2 to G4UrbanMscModel93, 
+//           now it is a frozen version of the Urban model corresponding 
+//           to g4 9.3
+
+// Class Description:
+//
+// Implementation of the model of multiple scattering based on
+// H.W.Lewis Phys Rev 78 (1950) 526 and others
+
+// -------------------------------------------------------------------
+// In its present form the model can be  used for simulation 
+//   of the e-/e+, muon and charged hadron multiple scattering
+//
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4UrbanMscModel93.hh"
+#include "Randomize.hh"
+#include "G4Electron.hh"
+#include "G4LossTableManager.hh"
+#include "G4ParticleChangeForMSC.hh"
+
+#include "G4Poisson.hh"
+#include "globals.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+using namespace std;
+
+G4UrbanMscModel93::G4UrbanMscModel93(const G4String& nam)
+  : G4VMscModel(nam),
+    isInitialized(false)
+{
+  masslimite    = 0.6*MeV;
+  lambdalimit   = 1.*mm;
+  fr            = 0.02;
+  //facsafety     = 0.3;
+  taubig        = 8.0;
+  tausmall      = 1.e-16;
+  taulim        = 1.e-6;
+  currentTau    = taulim;
+  tlimitminfix  = 1.e-6*mm;            
+  stepmin       = tlimitminfix;
+  smallstep     = 1.e10;
+  currentRange  = 0. ;
+  rangeinit     = 0.;
+  tlimit        = 1.e10*mm;
+  tlimitmin     = 10.*tlimitminfix;            
+  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(); 
+  inside        = false;  
+  insideskin    = false;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4UrbanMscModel93::~G4UrbanMscModel93()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4UrbanMscModel93::Initialise(const G4ParticleDefinition* p,
+				   const G4DataVector&)
+{
+  skindepth = skin*stepmin;
+  if(isInitialized) return;
+  // set values of some data members
+  SetParticle(p);
+
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
+
+  isInitialized = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::ComputeCrossSectionPerAtom( 
+                             const G4ParticleDefinition* part,
+                                   G4double KineticEnergy,
+                                   G4double AtomicNumber,G4double,
+				   G4double, G4double)
+{
+  const G4double sigmafactor = twopi*classic_electr_radius*classic_electr_radius;
+  const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
+                            Bohr_radius*Bohr_radius/(hbarc*hbarc);
+  const G4double epsmin = 1.e-4 , epsmax = 1.e10;
+
+  const G4double Zdat[15] = { 4.,  6., 13., 20., 26., 29., 32., 38., 47.,
+                             50., 56., 64., 74., 79., 82. };
+
+  const G4double Tdat[22] = { 100*eV,  200*eV,  400*eV,  700*eV,
+                               1*keV,   2*keV,   4*keV,   7*keV,
+			      10*keV,  20*keV,  40*keV,  70*keV,
+                             100*keV, 200*keV, 400*keV, 700*keV,
+                               1*MeV,   2*MeV,   4*MeV,   7*MeV,
+			      10*MeV,  20*MeV};
+
+  // corr. factors for e-/e+ lambda for T <= Tlim
+          G4double celectron[15][22] =
+          {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
+            1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
+            1.112,1.108,1.100,1.093,1.089,1.087            },
+           {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
+            1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
+            1.109,1.105,1.097,1.090,1.086,1.082            },
+           {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
+            1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
+            1.131,1.124,1.113,1.104,1.099,1.098            },
+           {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
+            1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
+            1.112,1.105,1.096,1.089,1.085,1.098            },
+           {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
+            1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
+            1.073,1.070,1.064,1.059,1.056,1.056            },
+           {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
+            1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
+            1.074,1.070,1.063,1.059,1.056,1.052            },
+           {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
+            1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
+            1.068,1.064,1.059,1.054,1.051,1.050            },
+           {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
+            1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
+            1.039,1.037,1.034,1.031,1.030,1.036            },
+           {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
+            1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
+            1.031,1.028,1.024,1.022,1.021,1.024            },
+           {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
+            1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
+            1.020,1.017,1.015,1.013,1.013,1.020            },
+           {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
+            1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
+            0.995,0.993,0.993,0.993,0.993,1.011            },
+           {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
+            1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
+            0.974,0.972,0.973,0.974,0.975,0.987            },
+           {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
+            1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
+            0.950,0.947,0.949,0.952,0.954,0.963            },
+           {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
+            1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
+            0.941,0.938,0.940,0.944,0.946,0.954            },
+           {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
+            1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
+            0.933,0.930,0.933,0.936,0.939,0.949            }};
+	    
+           G4double cpositron[15][22] = {
+           {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
+            1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
+            1.131,1.126,1.117,1.108,1.103,1.100            },
+           {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
+            1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
+            1.138,1.132,1.122,1.113,1.108,1.102            },
+           {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
+            1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
+            1.203,1.190,1.173,1.159,1.151,1.145            },
+           {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
+            1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
+            1.225,1.210,1.191,1.175,1.166,1.174            },
+           {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
+            1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
+            1.217,1.203,1.184,1.169,1.160,1.151            },
+           {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
+            1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
+            1.237,1.222,1.201,1.184,1.174,1.159            },
+           {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
+            1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
+            1.252,1.234,1.212,1.194,1.183,1.170            },
+           {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
+            2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
+            1.254,1.237,1.214,1.195,1.185,1.179            },
+           {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
+            2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
+            1.312,1.288,1.258,1.235,1.221,1.205            },
+           {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
+            2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
+            1.320,1.294,1.264,1.240,1.226,1.214            },
+           {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
+            2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
+            1.328,1.302,1.270,1.245,1.231,1.233            },
+           {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
+            2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
+            1.361,1.330,1.294,1.267,1.251,1.239            },
+           {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
+            3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
+            1.409,1.372,1.330,1.298,1.280,1.258            },
+           {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
+            3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
+            1.442,1.400,1.354,1.319,1.299,1.272            },
+           {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
+            3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
+            1.456,1.412,1.364,1.328,1.307,1.282            }};
+
+  //data/corrections for T > Tlim  
+  G4double Tlim = 10.*MeV;
+  G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
+                      ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
+  G4double bg2lim   = Tlim*(Tlim+2.*electron_mass_c2)/
+                      (electron_mass_c2*electron_mass_c2);
+
+  G4double sig0[15] = {0.2672*barn,  0.5922*barn, 2.653*barn,  6.235*barn,
+                      11.69*barn  , 13.24*barn  , 16.12*barn, 23.00*barn ,
+		      35.13*barn  , 39.95*barn  , 50.85*barn, 67.19*barn ,
+                      91.15*barn  , 104.4*barn  , 113.1*barn};
+		      		       
+  G4double hecorr[15] = {120.70, 117.50, 105.00, 92.92, 79.23,  74.510,  68.29,
+                          57.39,  41.97,  36.14, 24.53, 10.21,  -7.855, -16.84,
+			 -22.30};
+
+  G4double sigma;
+  SetParticle(part);
+
+  G4double Z23 = 2.*log(AtomicNumber)/3.; Z23 = exp(Z23);
+
+  // correction if particle .ne. e-/e+
+  // compute equivalent kinetic energy
+  // lambda depends on p*beta ....
+
+  G4double eKineticEnergy = KineticEnergy;
+
+  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 eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
+  G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
+                                 /(eTotalEnergy*eTotalEnergy);
+  G4double bg2   = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
+                                 /(electron_mass_c2*electron_mass_c2);
+
+  G4double eps = epsfactor*bg2/Z23;
+
+  if     (eps<epsmin)  sigma = 2.*eps*eps;
+  else if(eps<epsmax)  sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
+  else                 sigma = log(2.*eps)-1.+1./eps;
+
+  sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
+
+  // interpolate in AtomicNumber and beta2 
+  G4double c1,c2,cc1,cc2,corr;
+
+  // get bin number in Z
+  G4int iZ = 14;
+  while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
+  if (iZ==14)                               iZ = 13;
+  if (iZ==-1)                               iZ = 0 ;
+
+  G4double Z1 = Zdat[iZ];
+  G4double Z2 = Zdat[iZ+1];
+  G4double ratZ = (AtomicNumber-Z1)*(AtomicNumber+Z1)/
+                  ((Z2-Z1)*(Z2+Z1));
+
+  if(eKineticEnergy <= Tlim) 
+  {
+    // get bin number in T (beta2)
+    G4int iT = 21;
+    while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
+    if(iT==21)                                  iT = 20;
+    if(iT==-1)                                  iT = 0 ;
+
+    //  calculate betasquare values
+    G4double T = Tdat[iT],   E = T + electron_mass_c2;
+    G4double b2small = T*(E+electron_mass_c2)/(E*E);
+
+    T = Tdat[iT+1]; E = T + electron_mass_c2;
+    G4double b2big = T*(E+electron_mass_c2)/(E*E);
+    G4double ratb2 = (beta2-b2small)/(b2big-b2small);
+
+    if (charge < 0.)
+    {
+       c1 = celectron[iZ][iT];
+       c2 = celectron[iZ+1][iT];
+       cc1 = c1+ratZ*(c2-c1);
+
+       c1 = celectron[iZ][iT+1];
+       c2 = celectron[iZ+1][iT+1];
+       cc2 = c1+ratZ*(c2-c1);
+
+       corr = cc1+ratb2*(cc2-cc1);
+
+       sigma *= sigmafactor/corr;
+    }
+    else              
+    {
+       c1 = cpositron[iZ][iT];
+       c2 = cpositron[iZ+1][iT];
+       cc1 = c1+ratZ*(c2-c1);
+
+       c1 = cpositron[iZ][iT+1];
+       c2 = cpositron[iZ+1][iT+1];
+       cc2 = c1+ratZ*(c2-c1);
+
+       corr = cc1+ratb2*(cc2-cc1);
+
+       sigma *= sigmafactor/corr;
+    }
+  }
+  else
+  {
+    c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
+    c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
+    if((AtomicNumber >= Z1) && (AtomicNumber <= Z2))
+      sigma = c1+ratZ*(c2-c1) ;
+    else if(AtomicNumber < Z1)
+      sigma = AtomicNumber*AtomicNumber*c1/(Z1*Z1);
+    else if(AtomicNumber > Z2)
+      sigma = AtomicNumber*AtomicNumber*c2/(Z2*Z2);
+  }
+  return sigma;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::ComputeTruePathLengthLimit(
+                             const G4Track& track,
+			     G4PhysicsTable* theTable,
+			     G4double currentMinimalStep)
+{
+  tPathLength = currentMinimalStep;
+  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
+  G4StepStatus stepStatus = sp->GetStepStatus();
+
+  const G4DynamicParticle* dp = track.GetDynamicParticle();
+
+  if(stepStatus == fUndefined) {
+    inside = false;
+    insideskin = false;
+    tlimit = geombig;
+    SetParticle( dp->GetDefinition() );
+  }
+
+  theLambdaTable = theTable;
+  couple = track.GetMaterialCutsCouple();
+  currentMaterialIndex = couple->GetIndex();
+  currentKinEnergy = dp->GetKineticEnergy();
+  currentRange = 
+    theManager->GetRangeFromRestricteDEDX(particle,currentKinEnergy,couple);
+  lambda0 = GetLambda(currentKinEnergy);
+
+  // stop here if small range particle
+  if(inside) return tPathLength;            
+  
+  if(tPathLength > currentRange) tPathLength = currentRange;
+
+  presafety = sp->GetSafety();
+
+  //G4cout << "G4Urban2::StepLimit tPathLength= " 
+  //	 <<tPathLength<<" safety= " << presafety
+  //       << " range= " <<currentRange<< " lambda= "<<lambda0
+  //	 << " Alg: " << steppingAlgorithm <<G4endl;
+
+  // far from geometry boundary
+  if(currentRange < presafety)
+    {
+      inside = true;
+      return tPathLength;  
+    }
+
+  // standard  version
+  //
+  if (steppingAlgorithm == fUseDistanceToBoundary)
+    {
+      //compute geomlimit and presafety 
+      G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
+
+      // is it far from boundary ?
+      if(currentRange < presafety)
+	{
+	  inside = true;
+	  return tPathLength;   
+	}
+
+      smallstep += 1.;
+      insideskin = false;
+
+      if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
+        {
+          rangeinit = currentRange;
+          if(stepStatus == fUndefined) smallstep = 1.e10;
+          else  smallstep = 1.;
+
+          //define stepmin here (it depends on lambda!)
+          //rough estimation of lambda_elastic/lambda_transport
+          G4double rat = currentKinEnergy/MeV ;
+          rat = 1.e-3/(rat*(10.+rat)) ;
+          //stepmin ~ lambda_elastic
+          stepmin = rat*lambda0;
+          skindepth = skin*stepmin;
+          //define tlimitmin
+          tlimitmin = 10.*stepmin;
+          if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+	  //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
+	  //	 << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
+          // constraint from the geometry
+          if((geomlimit < geombig) && (geomlimit > geommin))
+            {
+              // 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 
+      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  
+      //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
+
+      // shortcut
+      if((tPathLength < tlimit) && (tPathLength < presafety) &&
+         (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
+	return tPathLength;   
+
+      // step reduction near to boundary
+      if(smallstep < skin)
+	{
+	  tlimit = stepmin;
+	  insideskin = true;
+	}
+      else if(geomlimit < geombig)
+	{
+	  if(geomlimit > skindepth)
+	    {
+	      if(tlimit > geomlimit-0.999*skindepth)
+		tlimit = geomlimit-0.999*skindepth;
+	    }
+	  else
+	    {
+	      insideskin = true;
+	      if(tlimit > stepmin) tlimit = stepmin;
+	    }
+	}
+
+      if(tlimit < stepmin) tlimit = stepmin;
+
+      if(tPathLength > tlimit) tPathLength = tlimit  ; 
+
+    }
+    // for 'normal' simulation with or without magnetic field 
+    //  there no small step/single scattering at boundaries
+  else if(steppingAlgorithm == fUseSafety)
+    {
+      // compute presafety again if presafety <= 0 and no boundary
+      // i.e. when it is needed for optimization purposes
+      if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix)) 
+	presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
+
+      // is far from boundary
+      if(currentRange < presafety)
+        {
+          inside = true;
+          return tPathLength;  
+        }
+
+      if((stepStatus == fGeomBoundary) || (stepStatus == fUndefined))
+      {
+        rangeinit = currentRange;
+        fr = facrange;
+        // 9.1 like stepping for e+/e- only (not for muons,hadrons)
+        if(mass < masslimite) 
+        {
+          if(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;
+    }
+  
+  // version similar to 7.1 (needed for some experiments)
+  else
+    {
+      if (stepStatus == fGeomBoundary)
+	{
+	  if (currentRange > lambda0) tlimit = facrange*currentRange;
+	  else                        tlimit = facrange*lambda0;
+
+	  if(tlimit < tlimitmin) tlimit = tlimitmin;
+	  if(tPathLength > tlimit) tPathLength = tlimit;
+	}
+    }
+  //G4cout << "tPathLength= " << tPathLength 
+  //	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
+  return tPathLength ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::ComputeGeomPathLength(G4double)
+{
+  lambdaeff = lambda0;
+  par1 = -1. ;  
+  par2 = par3 = 0. ;  
+
+  //  do the true -> geom transformation
+  zPathLength = tPathLength;
+
+  // z = t for very small tPathLength
+  if(tPathLength < tlimitminfix) return zPathLength;
+
+  // this correction needed to run MSC with eIoni and eBrem inactivated
+  // and makes no harm for a normal run
+  if(tPathLength > currentRange)
+    tPathLength = currentRange ;
+
+  G4double tau   = tPathLength/lambda0 ;
+
+  if ((tau <= tausmall) || insideskin) {
+    zPathLength  = tPathLength;
+    if(zPathLength > lambda0) zPathLength = lambda0;
+    return zPathLength;
+  }
+
+  G4double zmean = tPathLength;
+  if (tPathLength < currentRange*dtrl) {
+    if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
+    else             zmean = lambda0*(1.-exp(-tau));
+  } else if(currentKinEnergy < mass)  {
+    par1 = 1./currentRange ;
+    par2 = 1./(par1*lambda0) ;
+    par3 = 1.+par2 ;
+    if(tPathLength < currentRange)
+      zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
+    else
+      zmean = 1./(par1*par3) ;
+  } else {
+    G4double T1 = theManager->GetEnergy(particle,currentRange-tPathLength,couple);
+    G4double lambda1 = GetLambda(T1);
+
+    par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
+    par2 = 1./(par1*lambda0) ;
+    par3 = 1.+par2 ;
+    zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
+  }
+
+  zPathLength = zmean ;
+
+  //  sample z
+  if(samplez)
+  {
+    const G4double  ztmax = 0.99 ;
+    G4double zt = zmean/tPathLength ;
+
+    if (tPathLength > stepmin && zt < ztmax)              
+    {
+      G4double u,cz1;
+      if(zt >= third)
+      {
+        G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
+        cz1 = 1.+cz ;
+        G4double u0 = cz/cz1 ;
+        G4double grej ;
+        do {
+            u = exp(log(G4UniformRand())/cz1) ;
+            grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
+           } while (grej < G4UniformRand()) ;
+      }
+      else
+      {
+        cz1 = 1./zt-1.;
+        u = 1.-exp(log(G4UniformRand())/cz1) ;
+      }
+      zPathLength = tPathLength*u ;
+    }
+  }
+
+  if(zPathLength > lambda0) zPathLength = lambda0;
+  //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
+  return zPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::ComputeTrueStepLength(G4double geomStepLength)
+{
+  // step defined other than transportation 
+  if(geomStepLength == zPathLength && tPathLength <= currentRange)
+    return tPathLength;
+
+  // t = z for very small step
+  zPathLength = geomStepLength;
+  tPathLength = geomStepLength;
+  if(geomStepLength < tlimitminfix) return tPathLength;
+  
+  // recalculation
+  if((geomStepLength > lambda0*tausmall) && !insideskin)
+  {
+    if(par1 <  0.)
+      tPathLength = -lambda0*log(1.-geomStepLength/lambda0) ;
+    else 
+    {
+      if(par1*par3*geomStepLength < 1.)
+        tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
+      else 
+        tPathLength = currentRange;
+    }  
+  }
+  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
+  //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
+
+  return tPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::ComputeTheta0(G4double trueStepLength,
+					  G4double KineticEnergy)
+{
+  // for all particles take the width of the central part
+  //  from a  parametrization similar to the Highland formula
+  // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
+  const G4double c_highland = 13.6*MeV ;
+  G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
+                         KineticEnergy*(KineticEnergy+2.*mass)/
+                      ((currentKinEnergy+mass)*(KineticEnergy+mass)));
+  y = trueStepLength/currentRadLength;
+  G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
+  y = log(y);
+  // correction factor from e- scattering data
+  G4double corr = coeffth1+coeffth2*y;                
+
+  theta0 *= corr ;                                               
+
+  return theta0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4UrbanMscModel93::SampleScattering(const G4DynamicParticle* dynParticle,
+					 G4double safety)
+{
+  G4double kineticEnergy = dynParticle->GetKineticEnergy();
+
+  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) ||
+     (tPathLength/tausmall < lambda0)) return;
+
+  G4double cth  = SampleCosineTheta(tPathLength,kineticEnergy);
+
+  // protection against 'bad' cth values
+  if(std::abs(cth) > 1.) return;
+
+  G4double sth  = sqrt((1.0 - cth)*(1.0 + cth));
+  G4double phi  = twopi*G4UniformRand();
+  G4double dirx = sth*cos(phi);
+  G4double diry = sth*sin(phi);
+
+  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
+  G4ThreeVector newDirection(dirx,diry,cth);
+  newDirection.rotateUz(oldDirection);
+  fParticleChange->ProposeMomentumDirection(newDirection);
+
+  if (latDisplasment && safety > tlimitminfix) {
+
+    G4double r = SampleDisplacement();
+    /*    
+    G4cout << "G4UrbanMscModel93::SampleSecondaries: e(MeV)= " << kineticEnergy
+	   << " sinTheta= " << sth << " r(mm)= " << r
+           << " trueStep(mm)= " << tPathLength
+           << " geomStep(mm)= " << zPathLength
+           << G4endl;
+    */
+    if(r > 0.)
+      {
+        G4double latcorr = LatCorrelation();
+        if(latcorr > r) latcorr = r;
+
+        // sample direction of lateral displacement
+        // compute it from the lateral correlation
+        G4double Phi = 0.;
+        if(std::abs(r*sth) < latcorr)
+          Phi  = twopi*G4UniformRand();
+        else
+        {
+          G4double psi = std::acos(latcorr/(r*sth));
+          if(G4UniformRand() < 0.5)
+            Phi = phi+psi;
+          else
+            Phi = phi-psi;
+        }
+
+        dirx = std::cos(Phi);
+        diry = std::sin(Phi);
+
+        G4ThreeVector latDirection(dirx,diry,0.0);
+        latDirection.rotateUz(oldDirection);
+
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::SampleCosineTheta(G4double trueStepLength,
+					      G4double KineticEnergy)
+{
+  G4double cth = 1. ;
+  G4double tau = trueStepLength/lambda0 ;
+
+  Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
+         couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
+
+  if(Zold != Zeff)  
+    UpdateCache();
+
+  if(insideskin)
+  {
+    //no scattering, single or plural scattering
+    G4double mean = trueStepLength/stepmin ;
+
+    G4int n = G4Poisson(mean);
+    if(n > 0)
+    {
+      //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;
+      G4double sx=0.,sy=0.,sz=0.;
+      for(G4int i=1; i<=n; i++)
+      {
+        ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
+        if(ct < -1.) ct = -1.;
+        if(ct >  1.) ct =  1.; 
+        st = sqrt(1.-ct*ct);
+        phi = twopi*G4UniformRand();
+        sx += st*cos(phi);
+        sy += st*sin(phi);
+        sz += ct;
+      }
+      cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
+    }
+  }
+  else
+  {
+    if(trueStepLength >= currentRange*dtrl) 
+    {
+      if(par1*trueStepLength < 1.)
+	tau = -par2*log(1.-par1*trueStepLength) ;
+      // for the case if ioni/brems are inactivated
+      // see the corresponding condition in ComputeGeomPathLength 
+      else if(1.-KineticEnergy/currentKinEnergy > taulim)
+	tau = taubig ;
+    }
+    currentTau = tau ;
+    lambdaeff = trueStepLength/currentTau;
+    currentRadLength = couple->GetMaterial()->GetRadlen();
+
+    if (tau >= taubig) cth = -1.+2.*G4UniformRand();
+    else if (tau >= tausmall)
+    {
+      G4double xsi = 3.;
+      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);
+
+      // 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(-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 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.;
+
+      //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)
+      {
+        if(G4UniformRand() < prob)
+          cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
+        else
+          cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
+      }
+      else
+        cth = -1.+2.*G4UniformRand();
+    }
+  }  
+  return cth ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::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;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::SampleDisplacement()
+{
+  const G4double kappa = 2.5;
+  const G4double kappapl1 = kappa+1.;
+  const G4double kappami1 = kappa-1.;
+  G4double rmean = 0.0;
+  if ((currentTau >= tausmall) && !insideskin) {
+    if (currentTau < taulim) {
+      rmean = kappa*currentTau*currentTau*currentTau*
+             (1.-kappapl1*currentTau*0.25)/6. ;
+
+    } else {
+      G4double etau = 0.0;
+      if (currentTau<taubig) etau = exp(-currentTau);
+      rmean = -kappa*currentTau;
+      rmean = -exp(rmean)/(kappa*kappami1);
+      rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
+    }
+    if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean/3.0);
+    else          rmean = 0.;
+  }
+
+  // protection against z > t ...........................
+  if(rmean > 0.) {
+    G4double zt = (tPathLength-zPathLength)*(tPathLength+zPathLength);
+    if(zt <= 0.)
+      rmean = 0.;
+    else if(rmean*rmean > zt)
+      rmean = sqrt(zt);
+  }
+  return rmean;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4UrbanMscModel93::LatCorrelation()
+{
+  const G4double kappa = 2.5;
+  const G4double kappami1 = kappa-1.;
+
+  G4double latcorr = 0.;
+  if((currentTau >= tausmall) && !insideskin)
+  {
+    if(currentTau < taulim)
+      latcorr = lambdaeff*kappa*currentTau*currentTau*
+                (1.-(kappa+1.)*currentTau/3.)/3.;
+    else
+    {
+      G4double etau = 0.;
+      if(currentTau < taubig) etau = exp(-currentTau);
+      latcorr = -kappa*currentTau;
+      latcorr = exp(latcorr)/kappami1;
+      latcorr += 1.-kappa*etau/kappami1 ;
+      latcorr *= 2.*lambdaeff/3. ;
+    }
+  }
+
+  return latcorr;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4alphaIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4alphaIonisation.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/standard/src/G4alphaIonisation.cc	(revision 1197)
@@ -0,0 +1,158 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+// $Id: G4alphaIonisation.cc,v 1.1 2009/11/10 11:50:30 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4alphaIonisation
+//
+// Author:        Vladimir Ivanchenko
+//
+// Creation date: 28.10.2009 created from G4ionIonisation
+//
+// Modifications:
+//
+//
+//
+// -------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#include "G4alphaIonisation.hh"
+#include "G4Electron.hh"
+#include "G4Alpha.hh"
+#include "G4BraggIonModel.hh"
+#include "G4BetheBlochModel.hh"
+#include "G4UnitsTable.hh"
+#include "G4LossTableManager.hh"
+#include "G4IonFluctuations.hh"
+#include "G4UniversalFluctuation.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+G4alphaIonisation::G4alphaIonisation(const G4String& name)
+  : G4VEnergyLossProcess(name),
+    theParticle(0),
+    isInitialised(false),
+    nuclearStopping(true)
+{
+  //  SetLinearLossLimit(0.15);
+  SetStepFunction(0.2, 0.1*mm);
+  //  SetIntegral(true);
+  SetProcessSubType(fIonisation);
+  //  SetVerboseLevel(1);
+  mass = 0.0;
+  ratio = 0.0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4alphaIonisation::~G4alphaIonisation()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4alphaIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() == 2*eplus);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4alphaIonisation::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....
+
+void G4alphaIonisation::InitialiseEnergyLossProcess(
+		      const G4ParticleDefinition* part,
+		      const G4ParticleDefinition* bpart)
+{
+  if(!isInitialised) {
+
+    theParticle = part;
+    G4String pname = part->GetParticleName();
+
+    // define base particle
+    const G4ParticleDefinition* theBaseParticle = 0;
+    if(bpart == 0) { 
+      if(pname != "alpha") { theBaseParticle = G4Alpha::Alpha(); }
+    } else { theBaseParticle = bpart; }
+
+    mass  = part->GetPDGMass();
+    ratio = electron_mass_c2/mass;
+
+    SetBaseParticle(theBaseParticle);
+    SetSecondaryParticle(G4Electron::Electron());
+
+    if (!EmModel(1)) SetEmModel(new G4BraggIonModel(), 1);
+    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
+
+    // model limit defined for alpha
+    eth = (EmModel(1)->HighEnergyLimit())*mass/proton_mass_c2;
+    EmModel(1)->SetHighEnergyLimit(eth);
+
+    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
+    AddEmModel(1, EmModel(1), new G4IonFluctuations());
+
+    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);  
+    EmModel(2)->SetLowEnergyLimit(eth);
+    EmModel(2)->SetHighEnergyLimit(MaxKinEnergy());
+    AddEmModel(2, EmModel(2), FluctModel());    
+
+    isInitialised = true;
+  }
+  // reinitialisation of corrections for the new run
+  EmModel(1)->ActivateNuclearStopping(nuclearStopping);
+  EmModel(2)->ActivateNuclearStopping(nuclearStopping);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4alphaIonisation::PrintInfo()
+{
+  if (G4Alpha::Alpha() == theParticle) {
+    if(EmModel(1) && EmModel(2)) {
+      G4cout << "      NuclearStopping= " << nuclearStopping
+	     << G4endl;
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
