Index: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc	(revision 1058)
@@ -0,0 +1,749 @@
+//
+// ********************************************************************
+// * 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: G4GoudsmitSaundersonMscModel.cc,v 1.7 2009/04/20 19:22:29 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+// File name:     G4GoudsmitSaundersonMscModel
+//
+// Author:        Omrane Kadri
+//
+// Creation date: 20.02.2009
+//
+// Modifications:
+// 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
+//
+// 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta 
+//                     sampling from SampleCosineTheta() which means the splitting 
+//                     step into two sub-steps occur only for msc regime
+//
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//REFERENCES:
+//Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
+//Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
+//Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
+//Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
+//Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
+//Ref.6:G4UrbanMscModel G4_v9.1Ref09; 
+//Ref.7:G4eCoulombScatteringModel G4_v9.1Ref09.
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+#include "G4GoudsmitSaundersonMscModel.hh"
+#include "G4GoudsmitSaundersonTable.hh"
+
+#include "G4ParticleChangeForMSC.hh"
+#include "G4MaterialCutsCouple.hh"
+#include "G4DynamicParticle.hh"
+#include "G4DataInterpolation.hh" 
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+
+#include "G4LossTableManager.hh"
+#include "G4Track.hh"
+#include "G4PhysicsTable.hh"
+#include "Randomize.hh"
+#include "G4Poisson.hh"
+
+using namespace std;
+
+G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
+G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
+G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
+G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
+G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
+  : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(GeV),isInitialized(false)
+{ 
+  fr=0.02,rangeinit=0.,masslimite=0.6*MeV;
+  particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
+  tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
+  tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
+  theManager=G4LossTableManager::Instance();
+  inside=false;insideskin=false;
+  samplez=false;
+
+  GSTable = new G4GoudsmitSaundersonTable();
+
+  if(ener[0] < 0.0){ 
+    G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
+    LoadELSEPAXSections();
+  }
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
+{
+  delete GSTable;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
+					      const G4DataVector&)
+{ 
+  skindepth=skin*stepmin;
+  SetParticle(p);
+  if(isInitialized) return;
+  fParticleChange = GetParticleChangeForMSC();
+  InitialiseSafetyHelper();
+  isInitialized=true;
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
+                          G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
+{  
+  //Build cross section table : Taken from Ref.7
+  G4double cs=0.0;
+  G4double kinEnergy = kineticEnergy;
+  if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
+  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
+
+  //value0=Lambda0;value1=Lambda1
+  G4double value0,value1;
+  CalculateIntegrals(p,Z,kinEnergy,value0,value1);
+  
+  if(value1 > 0.0) cs = 1./value1;
+
+  return cs;
+}  
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4GoudsmitSaundersonMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
+				       G4double safety)
+{    
+  G4double kineticEnergy = dynParticle->GetKineticEnergy();
+  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)) return ;
+
+  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
+  G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
+  G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2,nu_interm;
+
+  ///////////////////////////////////////////
+  // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4
+  G4double  eloss = theManager->GetEnergy(particle,tPathLength,currentCouple);
+  if(eloss>0.5*kineticEnergy)return;
+  G4double ee       = kineticEnergy - 0.5*eloss;
+  G4double ttau     = ee/electron_mass_c2;
+  G4double ttau2    = ttau*ttau;
+  G4double epsilonpp= eloss/ee;
+  G4double temp2  = 0.166666*(4+ttau*(6+ttau*(7+ttau*(4+ttau))))*(epsilonpp/(ttau+1)/(ttau+2))*(epsilonpp/(ttau+1)/(ttau+2));
+  G4double cst1=epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
+
+  kineticEnergy *= (1 - cst1);
+  tPathLength *= (1 - temp2);
+  ///////////////////////////////////////////
+  // additivity rule for mixture xsection calculation
+  const G4Material* mat = currentCouple->GetMaterial();
+  G4int nelm = mat->GetNumberOfElements();
+  const G4ElementVector* theElementVector = mat->GetElementVector();
+  const G4double* theFraction = mat->GetFractionVector();
+  G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume();
+  G4double scrA,llambda0,llambda1;
+  scrA=0.0;llambda0 =0.;llambda1=0.;
+  for(G4int i=0;i<nelm;i++)
+    {
+      G4double l0,l1;
+      CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,l0,l1);
+      llambda0 += (theFraction[i]/l0);
+      llambda1 += (theFraction[i]/l1);
+    } 
+  if(llambda0>DBL_MIN)llambda0 =1./llambda0;
+  if(llambda1>DBL_MIN)llambda1 =1./llambda1;
+  G4double g1=llambda0/llambda1;
+  G4double x1,x0;
+
+  x0=g1/2.;
+  do
+    {  
+      x1 = x0-(x0*((1.+x0)*std::log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*std::log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)
+      delta = std::abs( x1 - x0 );    
+      x0 = x1;  // new approximation becomes the old approximation for the next iteration
+    } while (delta > 1e-10);
+  scrA = x1;
+
+  G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
+  G4double lambdan=0.;
+  if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0;
+  if((lambdan<=1.0e-12)||(lambdan>1.0e+5))return;
+  G4bool noscatt=false;
+  G4bool singlescatt=false;
+  G4bool mscatt=false;
+
+  G4double epsilon1=G4UniformRand();
+  if(epsilon1<(exp(-lambdan)))noscatt=true;// no scattering     
+  else if(epsilon1<((1.+lambdan)*exp(-lambdan)))//single scattering
+    {singlescatt=true;
+     ws=G4UniformRand();
+     ws= 1.-2.*scrA*ws/(1.-ws + scrA);
+     G4double phi0=twopi*G4UniformRand(); 
+     us=sqrt(1.-ws*ws)*cos(phi0);
+     vs=sqrt(1.-ws*ws)*sin(phi0);
+     G4double rr=G4UniformRand();
+     x_coord=(rr*us);
+     y_coord=(rr*vs);
+     z_coord=((1.-rr)+rr*ws);     
+    }
+  else
+  {mscatt=true;
+  // Ref.2 subsection 4.4 "The best solution found"
+  // Sample first substep scattering angle
+  SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
+  phi1  = twopi*G4UniformRand();
+  cosPhi1 = cos(phi1);
+  sinPhi1 = sin(phi1);
+
+  // Sample second substep scattering angle
+  SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
+  phi2  = twopi*G4UniformRand();
+  cosPhi2 = cos(phi2);
+  sinPhi2 = sin(phi2);
+
+  // Scattering direction
+  us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
+  vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
+  ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2; 
+  }
+  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
+  G4ThreeVector newDirection(us,vs,ws);
+  newDirection.rotateUz(oldDirection);
+  fParticleChange->ProposeMomentumDirection(newDirection);
+  
+  if((safety > tlimitminfix)&&(latDisplasment))
+    {  
+      // Scattering coordinates
+      if(mscatt)
+      {
+      if(scrA<DBL_MIN)scrA=DBL_MIN;
+      q1       = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
+      if(q1<DBL_MIN)q1=DBL_MIN;
+      Gamma    = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
+      Eta      = atomPerVolume*tPathLength/llambda1;      
+      delta    = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
+
+      nu = G4UniformRand(); 
+      nu = std::sqrt(nu);
+      nu0 = (1.0 - nu)/2.;
+      nu1 = nu*delta;
+      nu2 = nu*(1.0-delta);
+      nu_interm = 1.0 - nu0 - nu1 - nu2; 
+      x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu_interm*us);
+      y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu_interm*vs);
+      z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu_interm*ws)  ;
+      }
+      G4double r=sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
+
+      G4double check= 1.- tPathLength/zPathLength;
+      if(check<=0.)    return;
+      else if(r>check) {r=check;x_coord /=r;y_coord /=r;z_coord /=r;}
+
+      r *=tPathLength;
+      
+      if(r > tlimitminfix) {
+
+	G4ThreeVector latDirection = G4ThreeVector(x_coord*tPathLength,y_coord*tPathLength,z_coord*tPathLength);
+	latDirection.rotateUz(oldDirection);
+ 
+	ComputeDisplacement(fParticleChange, latDirection, r, safety);
+      }     
+    }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
+						     G4double &cost, G4double &sint)
+{
+  G4double u,Qn1,r1,tet;
+  G4double xi=0.;
+  Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
+
+  if((lambdan<1.)||(Qn1<0.001))//plural scatt. or small angle scatt.
+	    {
+	      G4double xi1,lambdai=0.;
+	      G4int i=0;
+	      do {xi1=G4UniformRand();
+	      lambdai -=std::log(xi1);
+	      xi +=2.*scrA*xi1/(1.-xi1 + scrA);
+	      i++;
+	      }while((lambdai<lambdan)&&(i<30));
+	    
+	    }
+	  else {
+	    if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
+	    else{// procedure described by Benedito in Ref.1
+	         do{r1=G4UniformRand();
+	         u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
+ 	         xi = 2.*u;
+       	         tet=acos(1.-xi);
+                 }while(tet*r1*r1>sin(tet));
+	       } 
+	      }
+
+
+  if(xi<0.)xi=0.;
+  if(xi>2.)xi=2.;
+  cost=(1. - xi);
+  sint=sqrt(xi*(2.-xi));
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+// Cubic spline log-log interpolation of Lambda0 and Lambda1
+// Screening parameter calculated according to Eq. 37 of Ref.1 
+void G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z, 
+						      G4double kinEnergy,G4double &Lam0,
+						      G4double &Lam1)
+{ 
+  //////// BEGIN OF: LAMBDA CALCULATION ////////////////////////////////
+  G4double TCSEForThisAtom[106],FTCSEForThisAtom[106],TCSPForThisAtom[106],FTCSPForThisAtom[106];
+  G4double summ00=0.0;
+  G4double summ10=0.0;
+  G4double InterpolatedValue=0.0;
+  
+  
+  G4int  iZ = G4int(Z);
+  if(iZ > 103) iZ = 103;
+  for(G4int i=0;i<106;i++)
+    {
+      TCSEForThisAtom[i]=TCSE[iZ-1][i];FTCSEForThisAtom[i]=FTCSE[iZ-1][i];
+      TCSPForThisAtom[i]=TCSP[iZ-1][i];FTCSPForThisAtom[i]=FTCSP[iZ-1][i];
+    }
+
+  G4double kineticE = kinEnergy;
+  if(kineticE<lowKEnergy)kineticE=lowKEnergy;
+  if(kineticE>highKEnergy)kineticE=highKEnergy;
+  kineticE /= eV;
+    
+  if(p==G4Electron::Electron())        
+    {
+      MyValue= new G4DataInterpolation(ener,TCSEForThisAtom,106,0.0,0);  
+      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
+      delete  MyValue;
+      summ00 = std::exp(InterpolatedValue);  
+      MyValue= new G4DataInterpolation(ener,FTCSEForThisAtom,106,0.0,0);  
+      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
+      delete  MyValue;
+      summ10 = std::exp(InterpolatedValue);
+    }
+  if(p==G4Positron::Positron())        
+    {
+      MyValue= new G4DataInterpolation(ener,TCSPForThisAtom,106,0.0,0);  
+      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
+      delete  MyValue;
+      summ00 = std::exp(InterpolatedValue);  
+      MyValue= new G4DataInterpolation(ener,FTCSPForThisAtom,106,0.0,0);  
+      InterpolatedValue = MyValue ->CubicSplineInterpolation(std::log(kineticE));
+      delete  MyValue;
+      summ10 = std::exp(InterpolatedValue);
+    }
+
+  summ00 *=barn;
+  summ10 *=barn;
+
+  Lam0=1./((1.+1./Z)*summ00);
+  Lam1=1./((1.+1./Z)*summ10);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+//t->g->t step transformations taken from Ref.6 
+G4double G4GoudsmitSaundersonMscModel::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;
+  currentCouple = track.GetMaterialCutsCouple();
+  currentMaterialIndex = currentCouple->GetIndex();
+  currentKinEnergy = dp->GetKineticEnergy();
+  currentRange = 
+    theManager->GetRangeFromRestricteDEDX(particle,currentKinEnergy,currentCouple);
+
+  lambda1 = GetLambda(currentKinEnergy);
+
+  // stop here if small range particle
+  if(inside) return tPathLength;            
+  
+  if(tPathLength > currentRange) tPathLength = currentRange;
+
+  G4double presafety = sp->GetSafety();
+
+  // 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, tPathLength);
+   
+      // is 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*lambda1;
+          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;
+
+      // 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(lambda1 > currentRange)
+            rangeinit = lambda1;
+          if(lambda1 > lambdalimit)
+            fr *= 0.75+0.25*lambda1/lambdalimit;
+        }
+
+        //lower limit for tlimit
+        G4double rat = currentKinEnergy/MeV ;
+        rat = 1.e-3/(rat*(10.+rat)) ;
+        tlimitmin = 10.*lambda1*rat;
+        if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
+      }
+      //step limit
+      tlimit = fr*rangeinit;               
+
+      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 > lambda1) tlimit = facrange*currentRange;
+	  else                        tlimit = facrange*lambda1;
+
+	  if(tlimit < tlimitmin) tlimit = tlimitmin;
+	  if(tPathLength > tlimit) tPathLength = tlimit;
+	}
+    }
+  return tPathLength ;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
+{
+  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/lambda1 ;
+
+  if ((tau <= tausmall) || insideskin) {
+    zPathLength  = tPathLength;
+    if(zPathLength > lambda1) zPathLength = lambda1;
+    return zPathLength;
+  }
+
+  G4double zmean = tPathLength;
+  if (tPathLength < currentRange*dtrl) {
+    if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
+    else             zmean = lambda1*(1.-exp(-tau));
+  } else if(currentKinEnergy < mass) {
+    par1 = 1./currentRange ;
+    par2 = 1./(par1*lambda1) ;
+    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,currentCouple);
+
+    lambda11 = GetLambda(T1);
+
+    par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
+    par2 = 1./(par1*lambda1) ;
+    par3 = 1.+par2 ;
+    zmean = (1.-exp(par3*log(lambda11/lambda1)))/(par1*par3) ;
+  }
+
+  zPathLength = zmean ;
+  //  sample z
+  if(samplez)
+  {
+    const G4double  ztmax = 0.99, onethird = 1./3. ;
+    G4double zt = zmean/tPathLength ;
+
+    if (tPathLength > stepmin && zt < ztmax)              
+    {
+      G4double u,cz1;
+      if(zt >= onethird)
+      {
+        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 > lambda1) zPathLength = lambda1;
+
+ return zPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4GoudsmitSaundersonMscModel::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 > lambda1*tausmall) && !insideskin)
+    {
+      if(par1 <  0.)
+	tPathLength = -lambda1*log(1.-geomStepLength/lambda1) ;
+      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......
+void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
+{ 
+  ///////////////////////////////////////
+  //Total & first transport x sections of e-/e+ from ELSEPA code
+  G4String filename = "XSECTIONS.dat";
+
+    char* path = getenv("G4LEDATA");
+    if (!path)
+      {
+        G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set";
+        G4Exception(excep);
+      }
+
+    G4String pathString(path);
+    G4String dirFile = pathString + "/msc_GS/" + filename;
+    FILE *infile;
+    infile = fopen(dirFile,"r"); 
+    if (infile == 0)
+      {
+	G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
+	G4Exception(excep);
+      }
+
+    // Read parameters from tables and take logarithms
+    G4float aRead;
+    for(G4int i=0 ; i<106 ;i++){
+	  fscanf(infile,"%f\t",&aRead);
+          if(aRead > 0.0) aRead = std::log(aRead);
+          else  aRead = 0.0;
+          ener[i]=aRead;
+    }        
+    for(G4int j=0;j<103;j++){
+      for(G4int i=0;i<106;i++){
+	  fscanf(infile,"%f\t",&aRead);
+          if(aRead > 0.0) aRead = std::log(aRead);
+          else  aRead = 0.0;
+	  TCSE[j][i]=aRead;
+	}        
+     }
+    for(G4int j=0;j<103;j++){
+      for(G4int i=0;i<106;i++){
+	  fscanf(infile,"%f\t",&aRead);
+          if(aRead > 0.0) aRead = std::log(aRead);
+          else  aRead = 0.0;
+	  FTCSE[j][i]=aRead;      
+	}        
+      }    
+    for(G4int j=0;j<103;j++){
+      for(G4int i=0;i<106;i++){
+	  fscanf(infile,"%f\t",&aRead);
+          if(aRead > 0.0) aRead = std::log(aRead);
+          else  aRead = 0.0;
+	  TCSP[j][i]=aRead;      
+	}        
+     }
+    for(G4int j=0;j<103;j++){
+      for(G4int i=0;i<106;i++){
+	  fscanf(infile,"%f\t",&aRead);
+          if(aRead > 0.0) aRead = std::log(aRead);
+          else  aRead = 0.0;
+	  FTCSP[j][i]=aRead;      
+	}        
+     }
+
+    fclose(infile);
+   //End loading XSections and Energies
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc	(revision 1058)
@@ -0,0 +1,284 @@
+//
+// ********************************************************************
+// * 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: G4GoudsmitSaundersonTable.cc,v 1.1 2009/03/05 18:48:30 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Header file
+//
+// File name:     G4GoudsmitSaundersonTable
+//
+// Author:        Omrane Kadri
+//
+// Creation date: 20.02.2009
+//
+// Modifications:
+// 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+  
+#include "G4GoudsmitSaundersonTable.hh"
+#include "G4ios.hh"
+
+using namespace std;
+
+static const G4double LAMBDAN[76]={1.0,1.1659144,1.3593564,1.5848932,1.8478498,2.1544347,2.5118864,2.9286446,
+  3.4145489,3.9810717,4.6415888,5.4116953,6.3095734,7.3564225,8.576959,10,11.659144,13.593564,15.848932,
+  18.478498,21.544347,25.118864,29.286446,34.145489,39.810717,46.415888,54.116953,63.095734,73.564225,
+  85.76959,100,116.59144,135.93564,158.48932,184.78498,215.44347,251.18864,292.86446,341.45489,398.10717,
+  464.15888,541.16953,630.95734,735.64225,857.6959,1000,1165.9144,1359.3564,1584.8932,1847.8498,2154.4347,
+  2511.8864,2928.6446,3414.5489,3981.0717,4641.5888,5411.6953,6309.5734,7356.4225,8576.959,10000,
+  11659.144,13593.564,15848.932,18478.498,21544.347,25118.864,29286.446,34145.489,39810.717,46415.888,
+  54116.953,63095.734,73564.225,85769.59,100000.};
+
+static const G4double scrA[76*11]={5.6991615e-05,0.00013145766,0.0003032221,0.00069941638,0.0016132837,0.0037212229,0.0085834253,0.01979865,0.045667845,0.1053381,0.24297434,
+4.7936556e-05,0.00010863608,0.0002461962,0.00055794142,0.0012644331,0.0028655179,0.0064939715,0.014716944,0.033352231,0.075584397,0.17129292,
+4.0337064e-05,9.0113601e-05,0.00020131513,0.00044974099,0.0010047281,0.0022445775,0.0050144194,0.011202288,0.025026077,0.055908629,0.12490071,
+3.3955914e-05,7.4956228e-05,0.00016546267,0.00036525176,0.00080627763,0.0017798234,0.003928884,0.0086728433,0.019144931,0.042261615,0.093290708,
+2.8595186e-05,6.2482352e-05,0.00013652803,0.00029832267,0.00065185455,0.0014243448,0.0031122867,0.00680055,0.014859647,0.0324693,0.070947545,
+2.40896e-05,5.2174455e-05,0.00011300203,0.00024474542,0.00053008182,0.0011480776,0.0024865636,0.0053855232,0.011664234,0.025262979,0.054715818,
+2.0301081e-05,4.3629599e-05,9.3765545e-05,0.00020151406,0.00043307928,0.00093074232,0.0020002834,0.0042988628,0.0092388016,0.019855357,0.042671682,
+1.7114155e-05,3.6528558e-05,7.7966778e-05,0.00016641277,0.00035519244,0.00075812493,0.0016181465,0.0034537819,0.0073717733,0.015734358,0.033583511,
+1.4432191e-05,3.0615355e-05,6.494509e-05,0.00013776959,0.00029225395,0.00061996534,0.0013151474,0.002789854,0.0059181845,0.012554388,0.026631925,
+1.2174324e-05,2.5683001e-05,5.4180958e-05,0.00011430036,0.00024112847,0.00050868556,0.001073125,0.0022638687,0.004775866,0.010075185,0.02125465,
+1.0272776e-05,2.156287e-05,4.5261124e-05,9.500448e-05,0.00019941731,0.00041858301,0.00087861852,0.0018442471,0.003871131,0.0081256222,0.017055929,
+8.6707436e-06,1.8116946e-05,3.7854164e-05,7.9093776e-05,0.00016526123,0.00034530242,0.00072148662,0.0015074987,0.003149819,0.0065813388,0.013751273,
+7.3205871e-06,1.5231781e-05,3.1692422e-05,6.5941709e-05,0.00013720343,0.00028547608,0.00059398363,0.0012358884,0.002571485,0.0053504309,0.011132521,
+6.1823301e-06,1.2813803e-05,2.6558521e-05,5.5046502e-05,0.0001140921,0.00023647292,0.00049012546,0.0010158583,0.002105518,0.0043640006,0.0090450433,
+5.2224117e-06,1.0785627e-05,2.2275102e-05,4.6003832e-05,9.5009779e-05,0.0001962197,0.0004052443,0.000836934,0.0017284846,0.0035697665,0.0073724887,
+4.4126431e-06,9.0831013e-06,1.8696896e-05,3.8486186e-05,7.9220983e-05,0.00016307057,0.00033566877,0.00069094949,0.0014222687,0.0029276356,0.0060263225,
+3.7293369e-06,7.652937e-06,1.570452e-05,3.2227096e-05,6.6132919e-05,0.00013571074,0.00027849072,0.00057148817,0.0011727455,0.0024065799,0.0049385198,
+3.1525796e-06,6.4507944e-06,1.3199587e-05,2.7008936e-05,5.5265564e-05,0.00011308415,0.00023139229,0.00047347387,0.00096882012,0.0019823954,0.0040563687,
+2.6656237e-06,5.4397244e-06,1.1100817e-05,2.2653381e-05,4.6228636e-05,9.4338538e-05,0.00019251616,0.00039286674,0.00080172111,0.0016360681,0.0033387157,
+2.2543784e-06,4.5888959e-06,9.3409186e-06,1.901389e-05,3.8703689e-05,7.8783224e-05,0.00016036705,0.00032643485,0.00066447387,0.0013525686,0.0027532185,
+1.9069837e-06,3.8725517e-06,7.8640717e-06,1.5969735e-05,3.2430074e-05,6.585643e-05,0.00013373603,0.00027158055,0.00055150431,0.0011199513,0.0022743086,
+1.6134537e-06,3.2691521e-06,6.6238996e-06,1.3421231e-05,2.7193867e-05,5.5099745e-05,0.00011164215,0.00022620741,0.00045833755,0.00092867565,0.0018816666,
+1.3653772e-06,2.7606661e-06,5.581811e-06,1.1285904e-05,2.2819053e-05,4.6138009e-05,9.3286778e-05,0.00018861722,0.00038136653,0.00077108777,0.001559068,
+1.1556674e-06,2.331988e-06,4.7056515e-06,9.4953989e-06,1.9160492e-05,3.8663407e-05,7.8017777e-05,0.00015742983,0.00031767312,0.00064102346,0.0012935028,
+9.7835084e-07,1.9704515e-06,3.9685959e-06,7.9929669e-06,1.6098268e-05,3.2422782e-05,6.5301237e-05,0.00013152022,0.00026488885,0.00053350047,0.0010744988,
+8.2839142e-07,1.6654293e-06,3.3482418e-06,6.7314314e-06,1.3533123e-05,2.72075e-05,5.4698983e-05,0.0001099689,0.00022108561,0.00044447883,0.00089359698,
+7.0154172e-07,1.4079982e-06,2.8258605e-06,5.6715182e-06,1.1382769e-05,2.2845283e-05,4.5850614e-05,9.2022444e-05,0.00018468957,0.00037067303,0.00074394288,
+5.9421855e-07,1.1906618e-06,2.3857812e-06,4.7804943e-06,9.5788859e-06,1.9193634e-05,3.8459125e-05,7.7062235e-05,0.00015441298,0.00030940405,0.00061996646,
+5.0339814e-07,1.007118e-06,2.0148796e-06,4.0310468e-06,8.0646696e-06,1.6134493e-05,3.2279297e-05,6.4579222e-05,0.00012919971,0.00025848199,0.00051712917,
+4.2652819e-07,8.5206673e-07,1.7021564e-06,3.4003635e-06,6.7928372e-06,1.3569913e-05,2.7108341e-05,5.4153784e-05,0.00010818192,0.00021611283,0.00043172422,
+3.6145372e-07,7.2104791e-07,1.4383863e-06,2.8693727e-06,5.7239836e-06,1.141852e-05,2.2778296e-05,4.5439406e-05,9.0645045e-05,0.00018082376,0.00036071725,
+3.0635477e-07,6.1030666e-07,1.2158264e-06,2.4221165e-06,4.8252353e-06,9.6126241e-06,1.9149852e-05,3.8149503e-05,7.5999782e-05,0.00015140346,0.0003016194,
+2.5969396e-07,5.166805e-07,1.0279744e-06,2.0452317e-06,4.0691409e-06,8.0958593e-06,1.6107316e-05,3.2046706e-05,6.3759311e-05,0.0001268539,0.0002523853,
+2.2017227e-07,4.3750469e-07,8.693663e-07,1.7275192e-06,3.4327563e-06,6.8212359e-06,1.3554489e-05,2.6934148e-05,5.3520891e-05,0.00010635145,0.00021133115,
+1.8669185e-07,3.7053306e-07,7.3540835e-07,1.4595876e-06,2.8968884e-06,5.7495434e-06,1.1411296e-05,2.2648349e-05,4.4950876e-05,8.9215386e-05,0.00017706852,
+1.5832463e-07,3.138715e-07,6.2223622e-07,1.2335555e-06,2.4454686e-06,4.8480322e-06,9.6110072e-06,1.9053392e-05,3.7772498e-05,7.4882289e-05,0.00014845079,
+1.3429812e-07,2.6594404e-07,5.2663605e-07,1.0428718e-06,2.0651483e-06,4.0895127e-06,8.0982631e-06,1.6036597e-05,3.1756493e-05,6.2885839e-05,0.00012452977,
+1.1392341e-07,2.2535736e-07,4.4579019e-07,8.8183895e-07,1.7444079e-06,3.4506969e-06,6.8259889e-06,1.3502816e-05,2.6710567e-05,5.2837454e-05,0.0001045203,
+9.665235e-08,1.9099697e-07,3.7743359e-07,7.4585536e-07,1.4739022e-06,2.9126127e-06,5.7556822e-06,1.1373938e-05,2.2476306e-05,4.4415955e-05,8.7771409e-05,
+8.2010037e-08,1.6190202e-07,3.1962263e-07,6.3099043e-07,1.2456844e-06,2.4591968e-06,4.8548804e-06,9.5843749e-06,1.8921216e-05,3.7353758e-05,7.3742789e-05,
+6.959459e-08,1.3726093e-07,2.7071877e-07,5.3393676e-07,1.0530798e-06,2.0769819e-06,4.096417e-06,8.0793348e-06,1.5934816e-05,3.1428127e-05,6.1985476e-05,
+5.9065875e-08,1.1638807e-07,2.2934024e-07,4.5191014e-07,8.9047944e-07,1.754671e-06,3.4575425e-06,6.8130148e-06,1.3424903e-05,2.6453489e-05,5.2126045e-05,
+5.0135964e-08,9.8704051e-08,1.9432138e-07,3.8256584e-07,7.5316789e-07,1.4827823e-06,2.9191944e-06,5.7470982e-06,1.131447e-05,2.2275109e-05,4.3853619e-05,
+4.2561071e-08,8.3719116e-08,1.6467843e-07,3.2392823e-07,6.3717816e-07,1.2533517e-06,2.4653867e-06,4.8495019e-06,9.5391399e-06,1.8763822e-05,3.6909094e-05,
+3.6134756e-08,7.1019192e-08,1.3958101e-07,2.7433229e-07,5.3917224e-07,1.0596883e-06,2.0827096e-06,4.0933542e-06,8.0450718e-06,1.5811771e-05,3.107643e-05,
+3.0682173e-08,6.0254106e-08,1.1832791e-07,2.323741e-07,4.563397e-07,8.9616668e-07,1.7599054e-06,3.4561284e-06,6.7871961e-06,1.3328796e-05,2.6175288e-05,
+2.6055209e-08,5.1127655e-08,1.0032685e-07,1.9686953e-07,3.8631344e-07,7.5805572e-07,1.4875187e-06,2.9189305e-06,5.7277633e-06,1.1239484e-05,2.2055032e-05,
+2.2128376e-08,4.338923e-08,8.5077425e-08,1.6681947e-07,3.2709893e-07,6.4137425e-07,1.257604e-06,2.4659048e-06,4.8351362e-06,9.4807153e-06,1.8589748e-05,
+1.8795341e-08,3.6826758e-08,7.2156716e-08,1.4138067e-07,2.7701501e-07,5.4277093e-07,1.0634813e-06,2.0837381e-06,4.0827839e-06,7.9996258e-06,1.5674112e-05,
+1.5965985e-08,3.1260736e-08,6.1207222e-08,1.198412e-07,2.3464409e-07,4.5942338e-07,8.995319e-07,1.7612461e-06,3.4484467e-06,6.7519154e-06,1.321997e-05,
+1.3563922e-08,2.6539201e-08,5.192666e-08,1.0159982e-07,1.9879044e-07,3.8895383e-07,7.6102795e-07,1.489029e-06,2.9134376e-06,5.7004385e-06,1.1153491e-05,
+1.1524398e-08,2.2533483e-08,4.4059382e-08,8.6148651e-08,1.6844517e-07,3.2935832e-07,6.4398939e-07,1.2591828e-06,2.4620613e-06,4.8140317e-06,9.4128041e-06,
+9.7925058e-09,1.9134608e-08,3.7389126e-08,7.3058552e-08,1.427568e-07,2.7894754e-07,5.4506498e-07,1.0650599e-06,2.0811326e-06,4.0665442e-06,7.946049e-06,
+8.321689e-09,1.6250266e-08,3.1732878e-08,6.1966712e-08,1.2100615e-07,2.3629604e-07,4.6142958e-07,9.0106147e-07,1.7595573e-06,3.435994e-06,6.7096736e-06,
+7.0724613e-09,1.3802258e-08,2.6935788e-08,5.2566522e-08,1.0258617e-07,2.0020197e-07,3.9070405e-07,7.6247826e-07,1.488014e-06,2.9039329e-06,5.6671686e-06,
+6.0113304e-09,1.1724321e-08,2.2866768e-08,4.4598666e-08,8.6983916e-08,1.6965085e-07,3.3088199e-07,6.4534243e-07,1.2586568e-06,2.4548468e-06,4.7878605e-06,
+5.109884e-09,9.9602974e-09,1.9414829e-08,3.7843808e-08,7.3765977e-08,1.4378625e-07,2.802713e-07,5.4631094e-07,1.0648812e-06,2.0756897e-06,4.0459796e-06,
+4.3440169e-09,8.4625895e-09,1.6485991e-08,3.2116397e-08,6.2566028e-08,1.2188503e-07,2.3744451e-07,4.625662e-07,9.011263e-07,1.7554863e-06,3.4198669e-06,
+3.6932746e-09,7.1908446e-09,1.400065e-08,2.7259414e-08,5.3074365e-08,1.0333635e-07,2.0119696e-07,3.917326e-07,7.6270751e-07,1.4849996e-06,2.8913098e-06,
+3.1402991e-09,6.1108495e-09,1.1891377e-08,2.3139967e-08,4.5029104e-08,8.7624161e-08,1.705118e-07,3.318066e-07,6.456774e-07,1.2564527e-06,2.444988e-06,
+2.6703582e-09,5.1935918e-09,1.010104e-08,1.964556e-08,3.8208739e-08,7.431235e-08,1.4453043e-07,2.8109789e-07,5.4670856e-07,1.063296e-06,2.0680092e-06,
+2.2709466e-09,4.4144654e-09,8.5812256e-09,1.668094e-08,3.2425876e-08,6.3032266e-08,1.2252765e-07,2.3818001e-07,4.6299523e-07,9.0001078e-07,1.74952e-06,
+1.9314485e-09,3.7525989e-09,7.2909004e-09,1.4165444e-08,2.752195e-08,5.347222e-08,1.0389083e-07,2.0184884e-07,3.9217081e-07,7.6194615e-07,1.4803803e-06,
+1.6428508e-09,3.1902859e-09,6.1952822e-09,1.2030747e-08,2.3362756e-08,4.5368618e-08,8.8102257e-08,1.7108759e-07,3.3223854e-07,6.4518089e-07,1.2528901e-06,
+1.3975005e-09,2.7125035e-09,5.2648821e-09,1.0218967e-08,1.9834686e-08,3.8498488e-08,7.4724326e-08,1.4503751e-07,2.8151315e-07,5.4640797e-07,1.0605603e-06,
+1.1888988e-09,2.306504e-09,4.474696e-09,8.6810621e-09,1.6841555e-08,3.2673188e-08,6.3387093e-08,1.2297311e-07,2.3857201e-07,4.6283778e-07,8.97921e-07,
+1.0115263e-09,1.9614686e-09,3.8035184e-09,7.37547e-09,1.4301904e-08,2.7733076e-08,5.3777699e-08,1.0428129e-07,2.022137e-07,3.9211616e-07,7.6035939e-07,
+8.6069493e-10,1.6682144e-09,3.2333632e-09,6.2669628e-09,1.214674e-08,2.3543031e-08,4.5631526e-08,8.8443846e-08,1.7142346e-07,3.3225603e-07,6.4398462e-07,
+7.3242257e-10,1.4189468e-09,2.7489734e-09,5.3256784e-09,1.0317616e-08,1.9988663e-08,3.8724708e-08,7.5022678e-08,1.4534395e-07,2.8157972e-07,5.4551385e-07,
+6.2332613e-10,1.2070514e-09,2.337417e-09,4.5263342e-09,8.7651034e-09,1.6973346e-08,3.286835e-08,6.3648523e-08,1.2325336e-07,2.3867624e-07,4.6218901e-07,
+5.3053132e-10,1.0269023e-09,1.9876836e-09,3.8473824e-09,7.4470363e-09,1.4414566e-08,2.7900995e-08,5.4005476e-08,1.045336e-07,2.023364e-07,3.9164459e-07,
+4.5159583e-10,8.7373189e-10,1.6904661e-09,3.2706549e-09,6.3279493e-09,1.2243096e-08,2.3687516e-08,4.5829781e-08,8.8669866e-08,1.7155537e-07,3.3191936e-07,
+3.8444428e-10,7.4348804e-10,1.4378533e-09,2.7807065e-09,5.3776896e-09,1.0400071e-08,2.0113002e-08,3.8897122e-08,7.5224283e-08,1.4547844e-07,2.8134501e-07,
+3.2731284e-10,6.3272904e-10,1.2231296e-09,2.3644341e-09,4.5706918e-09,8.8356126e-09,1.7080139e-08,3.3017647e-08,6.3826472e-08,1.2338307e-07,2.3851203e-07,
+2.7870239e-10,5.3853141e-10,1.0405942e-09,2.0107207e-09,3.8852779e-09,7.5074498e-09,1.4506505e-08,2.8030648e-08,5.4163097e-08,1.0465834e-07,2.0222935e-07,
+2.3733868e-10,4.5841161e-10,8.8540646e-10,1.7101325e-09,3.3030629e-09,6.3797539e-09,1.2322278e-08,2.3800062e-08,4.5969012e-08,8.8787587e-08,1.7149021e-07
+};
+
+static const G4double uvalues[320]={0,7.61544e-09,3.04617e-08,6.85389e-08,1.21847e-07,1.90386e-07,2.74156e-07,
+  3.73156e-07,4.87388e-07,6.1685e-07,7.61543e-07,9.21467e-07,1.09662e-06,1.28701e-06,1.49262e-06,
+  1.71347e-06,1.94955e-06,2.20086e-06,2.4674e-06,2.74917e-06,3.04617e-06,3.3584e-06,3.68587e-06,
+  4.02856e-06,4.38648e-06,4.75964e-06,5.14803e-06,5.55164e-06,5.97049e-06,6.40457e-06,6.85388e-06,
+  7.31842e-06,7.79819e-06,8.29319e-06,8.80342e-06,9.32888e-06,9.86957e-06,1.04255e-05,1.09966e-05,
+  1.1583e-05,1.21846e-05,1.28015e-05,1.34336e-05,1.40809e-05,1.47434e-05,1.54212e-05,1.61142e-05,
+  1.68224e-05,1.75459e-05,1.82845e-05,1.90385e-05,7.61524e-05,0.000171338,0.000304586,0.000475889,
+  0.000685233,0.000932601,0.00121797,0.00154133,0.00190265,0.0023019,0.00273905,0.00321407,0.00372692,
+  0.00427757,0.00486597,0.00549207,0.00615583,0.0068572,0.00759612,0.00837255,0.00918641,0.0100376,
+  0.0109262,0.011852,0.012815,0.013815,0.0148521,0.0159262,0.0170371,0.0181848,0.0193692,0.0205901,
+  0.0218476,0.0231415,0.0244717,0.0258382,0.0272407,0.0286793,0.0301537,0.0316639,0.0332098,0.0347912,
+  0.0364081,0.0380602,0.0397476,0.04147,0.0432273,0.0450194,0.0468461,0.0487074,0.050603,0.0525328,
+  0.0544967,0.0564946,0.0585262,0.0605914,0.0626901,0.0648222,0.0669873,0.0691854,0.0714163,0.0736799,
+  0.075976,0.0783043,0.0806647,0.0830571,0.0854812,0.0879369,0.090424,0.0929422,0.0954915,0.0980716,
+  0.100682,0.103323,0.105995,0.108696,0.111427,0.114188,0.116978,0.119797,0.122645,0.125522,0.128428,
+  0.131361,0.134323,0.137313,0.14033,0.143375,0.146447,0.149545,0.152671,0.155823,0.159001,0.162205,
+  0.165435,0.16869,0.17197,0.175276,0.178606,0.181961,0.18534,0.188743,0.192169,0.195619,0.199092,
+  0.202589,0.206107,0.209649,0.213212,0.216797,0.220404,0.224032,0.22768,0.23135,0.23504,0.238751,
+  0.242481,0.246231,0.25,0.253788,0.257595,0.261421,0.265264,0.269126,0.273005,0.276901,0.280814,
+  0.284744,0.288691,0.292653,0.296632,0.300625,0.304634,0.308658,0.312697,0.316749,0.320816,0.324896,
+  0.32899,0.333097,0.337216,0.341348,0.345492,0.349647,0.353814,0.357992,0.362181,0.366381,0.37059,
+  0.37481,0.379039,0.383277,0.387524,0.39178,0.396044,0.400316,0.404596,0.408882,0.413176,0.417476,
+  0.421783,0.426095,0.430413,0.434737,0.439065,0.443398,0.447736,0.452077,0.456422,0.46077,0.465122,
+  0.469476,0.473832,0.47819,0.48255,0.486912,0.491274,0.495637,0.5,0.508726,0.51745,0.526168,0.534878,
+  0.543578,0.552264,0.560935,0.569587,0.578217,0.586824,0.595404,0.603956,0.612476,0.620961,0.62941,
+  0.637819,0.646186,0.654508,0.662784,0.67101,0.679184,0.687303,0.695366,0.703368,0.711309,0.719186,
+  0.726995,0.734736,0.742405,0.75,0.757519,0.76496,0.77232,0.779596,0.786788,0.793893,0.800908,0.807831,
+  0.81466,0.821394,0.82803,0.834565,0.840999,0.847329,0.853553,0.85967,0.865677,0.871572,0.877355,
+  0.883022,0.888573,0.894005,0.899318,0.904508,0.909576,0.914519,0.919335,0.924024,0.928584,0.933013,
+  0.93731,0.941474,0.945503,0.949397,0.953154,0.956773,0.960252,0.963592,0.96679,0.969846,0.972759,
+  0.975528,0.978152,0.980631,0.982963,0.985148,0.987185,0.989074,0.990814,0.992404,0.993844,0.995134,
+  0.996273,0.997261,0.998097,0.998782,0.999315,0.999695,0.999924,1.0};
+
+G4double G4GoudsmitSaundersonTable::PDF[]  = {-1.};
+G4double G4GoudsmitSaundersonTable::CPDF[] = {-1.};
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable()
+{
+  if(PDF[0] < 0.0){ 
+    G4cout << "### G4GoudsmitSaundersonTable loading PDF data" << G4endl;
+    LoadPDFandCPDFdata();
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+G4double G4GoudsmitSaundersonTable::SampleTheta(G4double lambda, G4double Chia2, G4double rndm)
+{ 
+  //Benedito's procedure
+  G4double A[11],ThisPDF[320],ThisCPDF[320];
+  G4double PDF1[320],PDF2[320],PDF3[320],PDF4[320],CPDF1[320],CPDF2[320],CPDF3[320],CPDF4[320];
+  G4double coeff,Ckj,CkjPlus1,CkPlus1j,CkPlus1jPlus1,aa,b,m,F;
+  G4int Ind0,Ind1,Ind2,Ind3,KIndex=0,JIndex=0,IIndex=0;
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Find Lambda and Chia2 Index
+  for(G4int k=0;k<75;k++){if((lambda>=LAMBDAN[k])&&(lambda<LAMBDAN[k+1]))KIndex=k;}  
+  for(G4int j=0;j<11;j++){A[j]=scrA[11*KIndex+j];}  
+  for(G4int j=0;j<10;j++){if((Chia2>=A[j])&&(Chia2<A[j+1]))JIndex=j;}
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Calculate some necessary coefficients for PDF and CPDF interpolation 
+  coeff=(log(LAMBDAN[KIndex+1]/LAMBDAN[KIndex]))*(log(A[JIndex+1]/A[JIndex]));
+  Ckj=(log(LAMBDAN[KIndex+1]/lambda))*(log(A[JIndex+1]/Chia2))/coeff;
+  CkjPlus1=(log(LAMBDAN[KIndex+1]/lambda))*(log(Chia2/A[JIndex]))/coeff;
+  CkPlus1j=(log(lambda/LAMBDAN[KIndex]))*(log(A[JIndex+1]/Chia2))/coeff;
+  CkPlus1jPlus1=(log(lambda/LAMBDAN[KIndex]))*(log(Chia2/A[JIndex]))/coeff;
+  ///////////////////////////////////////////////////////////////////////////
+  // Calculate Interpolated PDF and CPDF arrays
+  Ind0=320*(11*KIndex+JIndex);
+  Ind1=320*(11*KIndex+JIndex+1);
+  Ind2=320*(11*(KIndex+1)+JIndex);
+  Ind3=320*(11*(KIndex+1)+JIndex+1);
+  for(G4int i=0 ; i<320 ;i++){
+    PDF1[i]=PDF[Ind0+i];PDF2[i]=PDF[Ind1+i];    
+    PDF3[i]=PDF[Ind2+i];PDF4[i]=PDF[Ind3+i];    
+    CPDF1[i]=CPDF[Ind0+i];CPDF2[i]=CPDF[Ind1+i];    
+    CPDF3[i]=CPDF[Ind2+i];CPDF4[i]=CPDF[Ind3+i];  
+    ThisPDF[i]=Ckj*PDF1[i]+CkjPlus1*PDF2[i]+CkPlus1j*PDF3[i]+CkPlus1jPlus1*PDF4[i];
+    ThisCPDF[i]=Ckj*CPDF1[i]+CkjPlus1*CPDF2[i]+CkPlus1j*CPDF3[i]+CkPlus1jPlus1*CPDF4[i];  
+  }
+
+  ///////////////////////////////////////////////////////////////////////////
+  // Find u Index using secant method
+  for(G4int i=0;i<319;i++){if((rndm>=ThisCPDF[i])&&(rndm<ThisCPDF[i+1]))IIndex=i;}
+  //CPDF^-1(rndm)=x ==> CPDF(x)=rndm;
+  aa=uvalues[IIndex];
+  b=uvalues[IIndex+1];
+  do{
+    m=0.5*(aa+b);
+    F=(ThisCPDF[IIndex]+(m-uvalues[IIndex])*ThisPDF[IIndex]
+       +((m-uvalues[IIndex])*(m-uvalues[IIndex])*(ThisPDF[IIndex+1]-ThisPDF[IIndex]))
+       /(2.*(uvalues[IIndex+1]-uvalues[IIndex])))-rndm;
+    if(F>0.)b=m;
+    else aa=m;
+  } while(sqrt((b-aa)*(b-aa))>1.0e-10);
+
+  return m;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void G4GoudsmitSaundersonTable::LoadPDFandCPDFdata()
+{ 
+  ///////////////////////////////////////
+  //Probability and their cumulative loading data
+  G4String filename;
+
+  for(G4int level = 0; level < 2; level++){
+    if(level == 0) filename = "PDF.dat";
+    if(level == 1) filename = "CPDF.dat";
+
+    char* path = getenv("G4LEDATA");
+    if (!path)
+      {
+        G4String excep = "G4GoudsmitSaundersonTable: G4LEDATA environment variable not set";
+        G4Exception(excep);
+      }
+
+    G4String pathString(path);
+    G4String dirFile = pathString + "/msc_GS/" + filename;
+    FILE *infile = fopen(dirFile,"r"); 
+    if (infile == 0)
+      {
+	G4String excep = "G4GoudsmitSaunderson - data files: " + dirFile + " not found";
+	G4Exception(excep);
+      }
+
+    // Read parameters into tables. 
+    G4float aRead;
+    for(G4int k=0 ; k<76 ;k++){
+      for(G4int j=0 ; j<11 ;j++){
+        for(G4int i=0 ; i<320 ;i++){
+	  fscanf(infile,"%f\t",&aRead);
+	  G4int idx = 320*(11*k+j)+i;
+	  if(level == 0)       PDF[idx]  = aRead;
+	  else if (level == 1) CPDF[idx] = aRead;
+        }
+      }
+    }
+      
+    fclose(infile);
+  } //End loading PDF and CPDF parameters
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Index: trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4HeatedKleinNishinaCompton.cc	(revision 1058)
@@ -0,0 +1,279 @@
+//
+// ********************************************************************
+// * 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: G4HeatedKleinNishinaCompton.cc,v 1.5 2009/04/12 17:09:57 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4HeatedKleinNishinaCompton
+//
+// Author:        Vladimir Grichine on base of M. Maire and V. Ivanchenko code
+//
+// Creation date: 15.03.2009
+//
+// Modifications:
+// 
+//
+// Class Description:
+//
+// -------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#include <CLHEP/Random/RandGamma.h>
+#include "globals.hh"
+#include "G4RandomDirection.hh"
+#include "Randomize.hh"
+
+#include "G4HeatedKleinNishinaCompton.hh"
+#include "G4Electron.hh"
+#include "G4Gamma.hh"
+#include "Randomize.hh"
+#include "G4DataVector.hh"
+#include "G4ParticleChangeForGamma.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+G4HeatedKleinNishinaCompton::G4HeatedKleinNishinaCompton(const G4ParticleDefinition*,
+                                             const G4String& nam)
+  : G4VEmModel(nam)
+{
+  theGamma = G4Gamma::Gamma();
+  theElectron = G4Electron::Electron();
+  lowestGammaEnergy = 1.0*eV;
+  fTemperature = 1.0*keV;
+  fParticleChange = 0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4HeatedKleinNishinaCompton::~G4HeatedKleinNishinaCompton()
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4HeatedKleinNishinaCompton::Initialise(const G4ParticleDefinition*,
+                                       const G4DataVector&)
+{
+  if(!fParticleChange) fParticleChange = GetParticleChangeForGamma();
+}
+
+////////////////////////////////////////////////////////////////////////////
+//
+//
+
+G4double G4HeatedKleinNishinaCompton::ComputeCrossSectionPerAtom(
+                                       const G4ParticleDefinition*,
+                                             G4double GammaEnergy,
+                                             G4double Z, G4double,
+                                             G4double, G4double)
+{
+  G4double CrossSection = 0.0 ;
+  if ( Z < 0.9999 )                 return CrossSection;
+  if ( GammaEnergy < 0.01*keV      ) return CrossSection;
+  //  if ( GammaEnergy > (100.*GeV/Z) ) return CrossSection;
+
+  static const G4double a = 20.0 , b = 230.0 , c = 440.0;
+  
+  static const G4double
+    d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527   *barn, d4=-1.9798e+1*barn,
+    e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
+    f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
+     
+  G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
+           p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
+
+  G4double T0  = 15.0*keV; 
+  if (Z < 1.5) T0 = 40.0*keV; 
+
+  G4double X   = max(GammaEnergy, T0) / electron_mass_c2;
+  CrossSection = p1Z*std::log(1.+2.*X)/X
+               + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
+		
+  //  modification for low energy. (special case for Hydrogen)
+  if (GammaEnergy < T0) {
+    G4double dT0 = 1.*keV;
+    X = (T0+dT0) / electron_mass_c2 ;
+    G4double sigma = p1Z*log(1.+2*X)/X
+                    + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
+    G4double   c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);             
+    G4double   c2 = 0.150; 
+    if (Z > 1.5) c2 = 0.375-0.0556*log(Z);
+    G4double    y = log(GammaEnergy/T0);
+    CrossSection *= exp(-y*(c1+c2*y));          
+  }
+  //  G4cout << "e= " << GammaEnergy << " Z= " << Z << " cross= " << CrossSection << G4endl;
+  return CrossSection;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//
+//
+
+void G4HeatedKleinNishinaCompton::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
+					      const G4MaterialCutsCouple*,
+					      const G4DynamicParticle* aDynamicGamma,
+					      G4double,
+					      G4double)
+{
+  // The scattered gamma energy is sampled according to Klein - Nishina formula.
+  // The random number techniques of Butcher & Messel are used 
+  // (Nuc Phys 20(1960),15).
+  // Note : Effects due to binding of atomic electrons are negliged.
+
+  // We start to prepare a heated electron from Maxwell distribution. 
+  // Then we try to boost to the electron rest frame and make scattering.
+  // The final step is to recover new gamma 4momentum in the lab frame
+
+  G4double eMomentumC2   = CLHEP::RandGamma::shoot(1.5,1.);
+  eMomentumC2          *= 2*electron_mass_c2*fTemperature; // electron (pc)^2
+  G4ThreeVector eMomDir = G4RandomDirection();
+  eMomDir              *= std::sqrt(eMomentumC2);
+  G4double eEnergy      = std::sqrt(eMomentumC2+electron_mass_c2*electron_mass_c2);
+  G4LorentzVector electron4v = G4LorentzVector(eMomDir,eEnergy);
+  G4ThreeVector bst = electron4v.boostVector();
+
+  G4LorentzVector gamma4v = aDynamicGamma->Get4Momentum();
+  gamma4v.boost(-bst);
+
+  G4ThreeVector gammaMomV = gamma4v.vect();
+  G4double  gamEnergy0    = gammaMomV.mag();
+ 
+ 
+  // G4double gamEnergy0 = aDynamicGamma->GetKineticEnergy();
+  G4double E0_m = gamEnergy0 / electron_mass_c2 ;
+
+  // G4ThreeVector gamDirection0 = /aDynamicGamma->GetMomentumDirection();
+
+  G4ThreeVector gamDirection0 = gammaMomV/gamEnergy0;
+  
+  // sample the energy rate of the scattered gamma in the electron rest frame
+  //
+
+  G4double epsilon, epsilonsq, onecost, sint2, greject ;
+
+  G4double epsilon0   = 1./(1. + 2.*E0_m);
+  G4double epsilon0sq = epsilon0*epsilon0;
+  G4double alpha1     = - log(epsilon0);
+  G4double alpha2     = 0.5*(1.- epsilon0sq);
+
+  do 
+  {
+    if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) 
+    {
+      epsilon   = exp(-alpha1*G4UniformRand());   // epsilon0**r
+      epsilonsq = epsilon*epsilon; 
+
+    } 
+    else 
+    {
+      epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
+      epsilon   = sqrt(epsilonsq);
+    };
+
+    onecost = (1.- epsilon)/(epsilon*E0_m);
+    sint2   = onecost*(2.-onecost);
+    greject = 1. - epsilon*sint2/(1.+ epsilonsq);
+
+  } while (greject < G4UniformRand());
+ 
+  //
+  // scattered gamma angles. ( Z - axis along the parent gamma)
+  //
+
+  G4double cosTeta = 1. - onecost; 
+  G4double sinTeta = sqrt (sint2);
+  G4double Phi     = twopi * G4UniformRand();
+  G4double dirx    = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
+
+  //
+  // update G4VParticleChange for the scattered gamma
+  //
+   
+  G4ThreeVector gamDirection1 ( dirx,diry,dirz );
+  gamDirection1.rotateUz(gamDirection0);
+  G4double gamEnergy1  = epsilon*gamEnergy0;
+  gamDirection1       *= gamEnergy1;
+
+  G4LorentzVector gamma4vfinal = G4LorentzVector(gamDirection1,gamEnergy1);
+
+  
+  // kinematic of the scattered electron
+  //
+
+  G4double eKinEnergy = gamEnergy0 - gamEnergy1;
+  G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
+  eDirection = eDirection.unit();
+  G4double eFinalMom = std::sqrt(eKinEnergy*(eKinEnergy+2*electron_mass_c2));
+  eDirection *= eFinalMom;
+  G4LorentzVector e4vfinal = G4LorentzVector(eDirection,gamEnergy1+electron_mass_c2);
+  
+  gamma4vfinal.boost(bst);
+  e4vfinal.boost(bst);
+
+  gamDirection1 = gamma4vfinal.vect();
+  gamEnergy1 = gamDirection1.mag(); 
+  gamDirection1 /= gamEnergy1;
+
+
+
+
+  fParticleChange->SetProposedKineticEnergy(gamEnergy1);
+
+  if( gamEnergy1 > lowestGammaEnergy ) 
+  {
+    gamDirection1 /= gamEnergy1;
+    fParticleChange->ProposeMomentumDirection(gamDirection1);
+  } 
+  else 
+  { 
+    fParticleChange->ProposeTrackStatus(fStopAndKill);
+    gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
+    fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
+  }
+
+  eKinEnergy = e4vfinal.t()-electron_mass_c2;
+
+  if( eKinEnergy > DBL_MIN ) 
+  {
+    // create G4DynamicParticle object for the electron.
+    eDirection = e4vfinal.vect();
+    G4double eFinMomMag = eDirection.mag();
+    eDirection /= eFinMomMag;
+    G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
+    fvect->push_back(dp);
+  }
+}
+
+//////////////////////////////////////////////////////////////////////////
+
+
Index: trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc	(revision 1058)
+++ trunk/source/processes/electromagnetic/standard/src/G4PairProductionRelModel.cc	(revision 1058)
@@ -0,0 +1,531 @@
+//
+// ********************************************************************
+// * 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: G4PairProductionRelModel.cc,v 1.3 2009/05/15 17:12:33 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
+//
+// -------------------------------------------------------------------
+//
+// GEANT4 Class file
+//
+//
+// File name:     G4PairProductionRelModel
+//
+// Author:        Andreas Schaelicke
+//
+// Creation date: 02.04.2009
+//
+// Modifications:
+//
+// Class Description:
+//
+// Main References:
+//  J.W.Motz et.al., Rev. Mod. Phys. 41 (1969) 581.
+//  S.Klein,  Rev. Mod. Phys. 71 (1999) 1501.
+//  T.Stanev et.al., Phys. Rev. D25 (1982) 1291.
+//  M.L.Ter-Mikaelian, High-energy Electromagnetic Processes in Condensed Media, Wiley, 1972.
+//
+// -------------------------------------------------------------------
+//
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+#include "G4PairProductionRelModel.hh"
+#include "G4Gamma.hh"
+#include "G4Electron.hh"
+#include "G4Positron.hh"
+
+#include "G4ParticleChangeForGamma.hh"
+#include "G4LossTableManager.hh"
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+
+const G4double G4PairProductionRelModel::facFel = log(184.15);
+const G4double G4PairProductionRelModel::facFinel = log(1194.); // 1440.
+
+const G4double G4PairProductionRelModel::preS1 = 1./(184.15*184.15);
+const G4double G4PairProductionRelModel::logTwo = log(2.);
+
+const G4double G4PairProductionRelModel::xgi[]={ 0.0199, 0.1017, 0.2372, 0.4083,
+						  0.5917, 0.7628, 0.8983, 0.9801 };
+const G4double G4PairProductionRelModel::wgi[]={ 0.0506, 0.1112, 0.1569, 0.1813,
+					    0.1813, 0.1569, 0.1112, 0.0506 };
+const G4double G4PairProductionRelModel::Fel_light[]  = {0., 5.31  , 4.79  , 4.74 ,  4.71} ;
+const G4double G4PairProductionRelModel::Finel_light[] = {0., 6.144 , 5.621 , 5.805 , 5.924} ;
+
+
+
+G4PairProductionRelModel::G4PairProductionRelModel(const G4ParticleDefinition*,
+					 const G4String& nam)
+  : G4VEmModel(nam),
+    theCrossSectionTable(0),
+    nbins(10),
+    fLPMconstant(fine_structure_const*electron_mass_c2*electron_mass_c2/(4.*pi*hbarc)*0.5),
+    fLPMflag(true),
+    lpmEnergy(0.),
+    use_completescreening(false)
+{
+  fParticleChange = 0;
+  theGamma    = G4Gamma::Gamma();
+  thePositron = G4Positron::Positron();
+  theElectron = G4Electron::Electron();
+
+  nist = G4NistManager::Instance();  
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4PairProductionRelModel::~G4PairProductionRelModel()
+{
+  if(theCrossSectionTable) {
+    theCrossSectionTable->clearAndDestroy();
+    delete theCrossSectionTable;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4PairProductionRelModel::Initialise(const G4ParticleDefinition*,
+					  const G4DataVector&)
+{
+  fParticleChange = GetParticleChangeForGamma();
+
+  if(theCrossSectionTable) {
+    theCrossSectionTable->clearAndDestroy();
+    delete theCrossSectionTable;
+  }
+
+  const G4ElementTable* theElementTable = G4Element::GetElementTable();
+  size_t nvect = G4Element::GetNumberOfElements();
+  theCrossSectionTable = new G4PhysicsTable(nvect);
+  G4PhysicsLogVector* ptrVector;
+  G4double emin = LowEnergyLimit();
+  G4double emax = HighEnergyLimit();
+  G4int n = nbins*G4int(log10(emax/emin));
+  G4bool spline = G4LossTableManager::Instance()->SplineFlag(); 
+  G4double e, value;
+
+  for(size_t j=0; j<nvect ; j++) { 
+
+    ptrVector  = new G4PhysicsLogVector(emin, emax, n);
+    ptrVector->SetSpline(spline);
+    G4double Z = (*theElementTable)[j]->GetZ();
+    G4VEmModel::SetCurrentElement((*theElementTable)[j]);
+    G4int   iz = G4int(Z);
+    indexZ[iz] = j;
+ 
+    for(G4int i=0; i<nbins; i++) {
+      e = ptrVector->GetLowEdgeEnergy( i ) ;
+      value = ComputeCrossSectionPerAtom(theGamma, e, Z);  
+      ptrVector->PutValue( i, value );
+    }
+
+    theCrossSectionTable->insert(ptrVector);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+/*
+G4double G4PairProductionRelModel::ComputeRelXSectionPerAtom(G4double k, G4double Z)
+{
+
+  G4double cross = 0.0;
+  
+
+
+
+}
+*/
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4PairProductionRelModel::ComputeXSectionPerAtom(G4double totalEnergy, G4double Z)
+{
+  G4double cross = 0.0;
+
+  // number of intervals and integration step 
+  G4double vcut = electron_mass_c2/totalEnergy ;
+
+  // limits by the screening variable
+  G4double dmax = DeltaMax();
+  G4double dmin = min(DeltaMin(totalEnergy),dmax);
+  G4double vcut1 = 0.5 - 0.5*sqrt(1. - dmin/dmax) ;
+  vcut = max(vcut, vcut1);
+
+
+  G4double vmax = 0.5;
+  G4int n = 1;  // needs optimisation
+
+  G4double delta = (vmax - vcut)*totalEnergy/G4double(n);
+
+  G4double e0 = vcut*totalEnergy;
+  G4double xs; 
+
+  // simple integration
+  for(G4int l=0; l<n; l++,e0 += delta) {
+    for(G4int i=0; i<8; i++) {
+
+      G4double eg = (e0 + xgi[i]*delta);
+      if (fLPMflag && totalEnergy>100.*GeV) 
+	xs = ComputeRelDXSectionPerAtom(eg,totalEnergy,Z);
+      else
+	xs = ComputeDXSectionPerAtom(eg,totalEnergy,Z);
+      cross += wgi[i]*xs;
+
+    }
+  }
+
+  cross *= delta*2.;
+
+  return cross;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double 
+G4PairProductionRelModel::ComputeDXSectionPerAtom(G4double eplusEnergy, 
+						  G4double totalEnergy, 
+						  G4double /*Z*/)
+{
+  // most simple case - complete screening:
+
+  //  dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
+  //     * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3)  + 1/9 * y*(1-y) ]
+  // y = E+/k
+  G4double yp=eplusEnergy/totalEnergy;
+  G4double ym=1.-yp;
+
+  G4double cross = 0.;
+  if (use_completescreening) 
+    cross = (yp*yp + ym*ym + 2./3.*ym*yp)*(Fel - fCoulomb) + 1./9.*yp*ym;
+  else {
+    G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
+    cross = (yp*yp + ym*ym)*(0.25*Phi1(delta) - lnZ/3. - fCoulomb)
+      + 2./3.*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
+  }
+  return cross/totalEnergy;
+
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double 
+G4PairProductionRelModel::ComputeRelDXSectionPerAtom(G4double eplusEnergy, 
+						     G4double totalEnergy, 
+						     G4double /*Z*/)
+{
+  // most simple case - complete screening:
+
+  //  dsig/dE+ = 4 * alpha * Z**2 * r0**2 / k
+  //     * [ (y**2 + (1-y**2) + 2/3*y*(1-y) ) * ( log (183 * Z**-1/3)  + 1/9 * y*(1-y) ]
+  // y = E+/k
+  G4double yp=eplusEnergy/totalEnergy;
+  G4double ym=1.-yp;
+
+  CalcLPMFunctions(totalEnergy,eplusEnergy); // gamma
+
+  G4double cross = 0.;
+  if (use_completescreening) 
+    cross = xiLPM*(2./3.*phiLPM*(yp*yp + ym*ym) + gLPM)*(Fel - fCoulomb);
+  else {
+    G4double delta = 0.25*DeltaMin(totalEnergy)/(yp*ym);
+    cross = (1./3.*gLPM + 2./3.*phiLPM)*(yp*yp + ym*ym)
+                             *(0.25*Phi1(delta) - lnZ/3. - fCoulomb) 
+           + 2./3.*gLPM*ym*yp*(0.25*Phi2(delta) - lnZ/3. - fCoulomb);
+    cross *= xiLPM;
+  }
+  return cross/totalEnergy;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void  G4PairProductionRelModel::CalcLPMFunctions(G4double k, G4double eplus)
+{
+  // *** calculate lpm variable s & sprime ***
+  // Klein eqs. (78) & (79)
+  G4double sprime = sqrt(0.125*k*lpmEnergy/(eplus*(k-eplus)));
+
+  G4double s1 = preS1*z23;
+  G4double logS1 = 2./3.*lnZ-2.*facFel;
+  G4double logTS1 = logTwo+logS1;
+
+  xiLPM = 2.;
+
+  if (sprime>1) 
+    xiLPM = 1.;
+  else if (sprime>sqrt(2.)*s1) {
+    G4double h  = log(sprime)/logTS1;
+    xiLPM = 1+h-0.08*(1-h)*(1-sqr(1-h))/logTS1;
+  }
+
+  G4double s = sprime/sqrt(xiLPM); 
+//   G4cout<<"k="<<k<<" y="<<eplus/k<<G4endl;
+//   G4cout<<"s="<<s<<G4endl;
+  
+  // *** calculate supression functions phi and G ***
+  // Klein eqs. (77)
+  G4double s2=s*s;
+  G4double s3=s*s2;
+  G4double s4=s2*s2;
+
+  if (s<0.1) {
+    // high suppression limit
+    phiLPM = 6.*s - 18.84955592153876*s2 + 39.47841760435743*s3 
+      - 57.69873135166053*s4;
+    gLPM = 37.69911184307752*s2 - 236.8705056261446*s3 + 807.7822389*s4;
+  }
+  else if (s<1.9516) {
+    // intermediate suppression
+    // using eq.77 approxim. valid s<2.      
+    phiLPM = 1.-exp(-6.*s*(1.+(3.-pi)*s)
+		+s3/(0.623+0.795*s+0.658*s2));
+    if (s<0.415827397755) {
+      // using eq.77 approxim. valid 0.07<s<2
+      G4double psiLPM = 1-exp(-4*s-8*s2/(1+3.936*s+4.97*s2-0.05*s3+7.50*s4));
+      gLPM = 3*psiLPM-2*phiLPM;
+    }
+    else {
+      // using alternative parametrisiation
+      G4double pre = -0.16072300849123999 + s*3.7550300067531581 + s2*-1.7981383069010097 
+	+ s3*0.67282686077812381 + s4*-0.1207722909879257;
+      gLPM = tanh(pre);
+    }
+  }
+  else {
+    // low suppression limit valid s>2.
+    phiLPM = 1. - 0.0119048/s4;
+    gLPM = 1. - 0.0230655/s4;
+  }
+
+  // *** make sure suppression is smaller than 1 ***
+  // *** caused by Migdal approximation in xi    ***
+  if (xiLPM*phiLPM>1. || s>0.57)  xiLPM=1./phiLPM;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double 
+G4PairProductionRelModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
+						     G4double gammaEnergy, G4double Z,
+						     G4double, G4double, G4double)
+{
+  //  static const G4double gammaEnergyLimit = 1.5*MeV;
+  G4double crossSection = 0.0 ;
+  if ( Z < 1. ) return crossSection;
+  if ( gammaEnergy <= 2.0*electron_mass_c2 ) return crossSection;
+
+  SetCurrentElement(Z);
+  // choose calculator according to parameters and switches
+  // in the moment only one calculator:
+  crossSection=ComputeXSectionPerAtom(gammaEnergy,Z);
+
+  G4double xi = Finel/(Fel - fCoulomb); // inelastic contribution
+  crossSection*=4.*Z*(Z+xi)*fine_structure_const*classic_electr_radius*classic_electr_radius;
+
+  return crossSection;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void 
+G4PairProductionRelModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
+					    const G4MaterialCutsCouple* couple,
+					    const G4DynamicParticle* aDynamicGamma,
+					    G4double,
+					    G4double)
+// The secondaries e+e- energies are sampled using the Bethe - Heitler
+// cross sections with Coulomb correction.
+// A modified version of the random number techniques of Butcher & Messel
+// is used (Nuc Phys 20(1960),15).
+//
+// GEANT4 internal units.
+//
+// Note 1 : Effects due to the breakdown of the Born approximation at
+//          low energy are ignored.
+// Note 2 : The differential cross section implicitly takes account of 
+//          pair creation in both nuclear and atomic electron fields.
+//          However triplet prodution is not generated.
+{
+  const G4Material* aMaterial = couple->GetMaterial();
+
+  G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
+  G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
+
+  G4double epsil ;
+  G4double epsil0 = electron_mass_c2/GammaEnergy ;
+  if(epsil0 > 1.0) return;
+
+  // do it fast if GammaEnergy < 2. MeV
+  static const G4double Egsmall=2.*MeV;
+
+  // select randomly one element constituing the material
+  const G4Element* anElement = SelectRandomAtom(aMaterial, theGamma, GammaEnergy);
+
+  if (GammaEnergy < Egsmall) {
+
+    epsil = epsil0 + (0.5-epsil0)*G4UniformRand();
+
+  } else {
+    // now comes the case with GammaEnergy >= 2. MeV
+
+    // Extract Coulomb factor for this Element
+    G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
+    if (GammaEnergy > 50.*MeV) FZ += 8.*(anElement->GetfCoulomb());
+
+    // limits of the screening variable
+    G4double screenfac = 136.*epsil0/(anElement->GetIonisation()->GetZ3());
+    G4double screenmax = exp ((42.24 - FZ)/8.368) - 0.952 ;
+    G4double screenmin = min(4.*screenfac,screenmax);
+
+    // limits of the energy sampling
+    G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ;
+    G4double epsilmin = max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin;
+
+    //
+    // sample the energy rate of the created electron (or positron)
+    //
+    //G4double epsil, screenvar, greject ;
+    G4double  screenvar, greject ;
+
+    G4double F10 = ScreenFunction1(screenmin) - FZ;
+    G4double F20 = ScreenFunction2(screenmin) - FZ;
+    G4double NormF1 = max(F10*epsilrange*epsilrange,0.); 
+    G4double NormF2 = max(1.5*F20,0.);
+
+    do {
+      if ( NormF1/(NormF1+NormF2) > G4UniformRand() ) {
+	epsil = 0.5 - epsilrange*pow(G4UniformRand(), 0.333333);
+	screenvar = screenfac/(epsil*(1-epsil));
+	if (fLPMflag && GammaEnergy>100.*GeV) {
+	  CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
+	  greject = xiLPM*((gLPM+2.*phiLPM)*Phi1(screenvar) - gLPM*Phi2(screenvar) - phiLPM*FZ)/F10;
+	}
+	else {
+	  greject = (ScreenFunction1(screenvar) - FZ)/F10;
+	}
+              
+      } else { 
+	epsil = epsilmin + epsilrange*G4UniformRand();
+	screenvar = screenfac/(epsil*(1-epsil));
+	if (fLPMflag && GammaEnergy>100.*GeV) {
+	  CalcLPMFunctions(GammaEnergy,GammaEnergy*epsil);
+	  greject = xiLPM*((0.5*gLPM+phiLPM)*Phi1(screenvar) + 0.5*gLPM*Phi2(screenvar) - 0.5*(gLPM+phiLPM)*FZ)/F20;
+	}
+	else {
+	  greject = (ScreenFunction2(screenvar) - FZ)/F20;
+	}
+      }
+
+    } while( greject < G4UniformRand() );
+
+  }   //  end of epsil sampling
+   
+  //
+  // fixe charges randomly
+  //
+
+  G4double ElectTotEnergy, PositTotEnergy;
+  if (G4UniformRand() > 0.5) {
+
+    ElectTotEnergy = (1.-epsil)*GammaEnergy;
+    PositTotEnergy = epsil*GammaEnergy;
+     
+  } else {
+    
+    PositTotEnergy = (1.-epsil)*GammaEnergy;
+    ElectTotEnergy = epsil*GammaEnergy;
+  }
+
+  //
+  // scattered electron (positron) angles. ( Z - axis along the parent photon)
+  //
+  //  universal distribution suggested by L. Urban 
+  // (Geant3 manual (1993) Phys211),
+  //  derived from Tsai distribution (Rev Mod Phys 49,421(1977))
+
+  G4double u;
+  const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
+
+  if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
+  else                            u= - log(G4UniformRand()*G4UniformRand())/a2;
+
+  G4double TetEl = u*electron_mass_c2/ElectTotEnergy;
+  G4double TetPo = u*electron_mass_c2/PositTotEnergy;
+  G4double Phi  = twopi * G4UniformRand();
+  G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
+  G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
+   
+  //
+  // kinematic of the created pair
+  //
+  // the electron and positron are assumed to have a symetric
+  // angular distribution with respect to the Z axis along the parent photon.
+
+  G4double ElectKineEnergy = max(0.,ElectTotEnergy - electron_mass_c2);
+
+  G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
+  ElectDirection.rotateUz(GammaDirection);   
+
+  // create G4DynamicParticle object for the particle1  
+  G4DynamicParticle* aParticle1= new G4DynamicParticle(
+		     theElectron,ElectDirection,ElectKineEnergy);
+  
+  // the e+ is always created (even with Ekine=0) for further annihilation.
+
+  G4double PositKineEnergy = max(0.,PositTotEnergy - electron_mass_c2);
+
+  G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
+  PositDirection.rotateUz(GammaDirection);   
+
+  // create G4DynamicParticle object for the particle2 
+  G4DynamicParticle* aParticle2= new G4DynamicParticle(
+                      thePositron,PositDirection,PositKineEnergy);
+
+  // Fill output vector
+  fvect->push_back(aParticle1);
+  fvect->push_back(aParticle2);
+
+  // kill incident photon
+  fParticleChange->SetProposedKineticEnergy(0.);
+  fParticleChange->ProposeTrackStatus(fStopAndKill);   
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+
+void G4PairProductionRelModel::SetupForMaterial(const G4ParticleDefinition*,
+						const G4Material* mat, G4double)
+{
+  lpmEnergy = mat->GetRadlen()*fLPMconstant;
+  //  G4cout<<" lpmEnergy="<<lpmEnergy<<G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
