Index: trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4EnergyLossForExtrapolator.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4EnergyLossForExtrapolator.cc,v 1.13 2007/07/28 13:44:25 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4EnergyLossForExtrapolator.cc,v 1.18 2008/11/13 14:14:07 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 //---------------------------------------------------------------------------
@@ -67,4 +67,6 @@
 #include "G4MuBremsstrahlungModel.hh"
 #include "G4ProductionCuts.hh"
+#include "G4LossTableManager.hh"
+#include "G4WentzelVIModel.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -88,4 +90,5 @@
   delete invRangePositron;
   delete invRangeProton;
+  delete mscElectron;
   delete cuts;
 }
@@ -100,6 +103,6 @@
   if(!isInitialised) Initialisation();
   G4double kinEnergyFinal = kinEnergy;
-  if(mat && part) {
-    G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
     G4double r  = ComputeRange(kinEnergy,part);
     if(r <= step) {
@@ -125,6 +128,6 @@
   G4double kinEnergyFinal = kinEnergy;
 
-  if(mat && part) {
-    G4double step = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    G4double step = TrueStepLength(kinEnergy,stepLength,mat,part);
     G4double r  = ComputeRange(kinEnergy,part);
 
@@ -141,8 +144,32 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4EnergyLossForExtrapolator::ComputeTrueStep(const G4Material* mat, 
-						      const G4ParticleDefinition* part,
-						      G4double kinEnergy, G4double stepLength)
-{
+G4double G4EnergyLossForExtrapolator::TrueStepLength(G4double kinEnergy, 
+						     G4double stepLength,
+						     const G4Material* mat, 
+						     const G4ParticleDefinition* part)
+{
+  G4double res = stepLength;
+  if(!isInitialised) Initialisation();
+  if(SetupKinematics(part, mat, kinEnergy)) {
+    if(part == electron || part == positron) {
+      G4double x = stepLength*ComputeValue(kinEnergy, mscElectron);
+      if(x < 0.2) res *= (1.0 + 0.5*x + x*x/3.0);
+      else if(x < 0.9999) res = -std::log(1.0 - x)*stepLength/x;
+      else res = ComputeRange(kinEnergy,part);
+    
+    } else {
+      res = ComputeTrueStep(mat,part,kinEnergy,stepLength);
+    }
+  }
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4EnergyLossForExtrapolator::SetupKinematics(const G4ParticleDefinition* part, 
+						    const G4Material* mat, 
+						    G4double kinEnergy)
+{
+  if(!part || !mat || kinEnergy < keV) return false;
   if(!isInitialised) Initialisation();
   G4bool flag = false;
@@ -182,6 +209,5 @@
     if(tmax > maxEnergyTransfer) tmax = maxEnergyTransfer;
   }
-  G4double theta = ComputeScatteringAngle(stepLength);
-  return stepLength*std::sqrt(1.0 + 0.625*theta*theta);
+  return true;
 } 
 
@@ -231,4 +257,5 @@
   invRangeMuon     = PrepareTable();
   invRangeProton   = PrepareTable();
+  mscElectron      = PrepareTable();
 
   G4LossTableBuilder builder; 
@@ -262,4 +289,5 @@
   builder.BuildInverseRangeTable(rangeProton, invRangeProton);  
 
+  ComputeTrasportXS(electron, mscElectron);
 }
 
@@ -273,4 +301,5 @@
 
     G4PhysicsVector* v = new G4PhysicsLogVector(emin, emax, nbins);
+    v->SetSpline(G4LossTableManager::Instance()->SplineFlag());
     table->push_back(v);
   }
@@ -488,2 +517,43 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+void G4EnergyLossForExtrapolator::ComputeTrasportXS(const G4ParticleDefinition* part, 
+						    G4PhysicsTable* table)
+{
+  G4DataVector v;
+  G4WentzelVIModel* msc = new G4WentzelVIModel();
+  msc->SetPolarAngleLimit(CLHEP::pi);
+  msc->Initialise(part, v);
+
+  mass    = part->GetPDGMass();
+  charge2 = 1.0;
+  currentParticle = part;
+
+  const G4MaterialTable* mtable = G4Material::GetMaterialTable();
+
+  if(0<verbose) {
+    G4cout << "G4EnergyLossForExtrapolator::ComputeProtonDEDX for " << part->GetParticleName() 
+           << G4endl;
+  }
+ 
+  for(G4int i=0; i<nmat; i++) {  
+
+    const G4Material* mat = (*mtable)[i];
+    if(1<verbose)
+      G4cout << "i= " << i << "  mat= " << mat->GetName() << G4endl;
+    G4PhysicsVector* aVector = (*table)[i];
+    for(G4int j=0; j<nbins; j++) {
+        
+       G4double e = aVector->GetLowEdgeEnergy(j);
+       G4double xs = msc->CrossSectionPerVolume(mat,part,e);
+       aVector->PutValue(j,xs);
+       if(1<verbose) {
+         G4cout << "j= " << j << "  e(MeV)= " << e/MeV  
+                << " xs(1/mm)= " << xs*mm << G4endl;
+       }
+    }
+  }
+  delete msc;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
Index: trunk/source/processes/electromagnetic/muons/src/G4MuBetheBlochModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuBetheBlochModel.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuBetheBlochModel.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuBetheBlochModel.cc,v 1.23 2007/05/22 17:35:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuBetheBlochModel.cc,v 1.25 2009/02/20 14:48:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -88,4 +88,5 @@
   theElectron = G4Electron::Electron();
   corr = G4LossTableManager::Instance()->EmCorrections();
+  fParticleChange = 0;
 
   if(p) SetParticle(p);
@@ -99,16 +100,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4MuBetheBlochModel::SetParticle(const G4ParticleDefinition* p)
-{
-  if(!particle) {
-    particle = p;
-    mass = particle->GetPDGMass();
-    massSquare = mass*mass;
-    ratio = electron_mass_c2/mass;
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
 G4double G4MuBetheBlochModel::MinEnergyCut(const G4ParticleDefinition*,
                                            const G4MaterialCutsCouple* couple)
@@ -117,4 +106,15 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4double G4MuBetheBlochModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
+						 G4double kinEnergy) 
+{
+  G4double tau  = kinEnergy/mass;
+  G4double tmax = 2.0*electron_mass_c2*tau*(tau + 2.) /
+                  (1. + 2.0*(tau + 1.)*ratio + ratio*ratio);
+  return tmax;
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -124,9 +124,12 @@
   if(p) SetParticle(p);
 
-  if(pParticleChange)
-    fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>
-                                                             (pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+  if(!fParticleChange) {
+    if(pParticleChange) {
+      fParticleChange = 
+	reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
+    } else {
+      fParticleChange = new G4ParticleChangeForLoss();
+    }
+  }
 }
 
@@ -275,5 +278,5 @@
 
   //High order corrections
-  dedx += corr->HighOrderCorrections(p,material,kineticEnergy);
+  dedx += corr->HighOrderCorrections(p,material,kineticEnergy,cutEnergy);
 
   return dedx;
Index: trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlung.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuBremsstrahlung.cc,v 1.38 2007/05/22 17:35:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuBremsstrahlung.cc,v 1.42 2009/02/20 14:48:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -82,5 +82,7 @@
     lowestKinEnergy(1.*GeV),
     isInitialised(false)
-{}
+{
+  SetProcessSubType(fBremsstrahlung);
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -91,6 +93,23 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4MuBremsstrahlung::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
-                                                     const G4ParticleDefinition*)
+G4bool G4MuBremsstrahlung::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4MuBremsstrahlung::MinPrimaryEnergy(const G4ParticleDefinition*,
+					      const G4Material*,
+					      G4double)
+{
+  return lowestKinEnergy;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4MuBremsstrahlung::InitialiseEnergyLossProcess(
+				 const G4ParticleDefinition* part,
+				 const G4ParticleDefinition*)
 {
   if(!isInitialised) {
@@ -105,7 +124,7 @@
     em->SetLowestKineticEnergy(lowestKinEnergy);
 
-    G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
-    em->SetLowEnergyLimit(0.1*keV);
-    em->SetHighEnergyLimit(100.0*TeV);
+    G4VEmFluctuationModel* fm = 0;
+    em->SetLowEnergyLimit(MinKinEnergy());
+    em->SetHighEnergyLimit(MaxKinEnergy());
     AddEmModel(1, em, fm);
   }
@@ -115,8 +134,5 @@
 
 void G4MuBremsstrahlung::PrintInfo()
-{
-  G4cout << "      Parametrised model "
-         << G4endl;
-}
+{}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuBremsstrahlungModel.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuBremsstrahlungModel.cc,v 1.24 2007/11/08 11:48:28 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuBremsstrahlungModel.cc,v 1.33 2009/02/20 14:48:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -46,9 +46,10 @@
 // 13-02-03 Add name (V.Ivanchenko)
 // 10-02-04 Add lowestKinEnergy (V.Ivanchenko)
-// 08-04-05 Major optimisation of internal interfaces (V.Ivantchenko)
-// 03-08-05 Angular correlations according to PRM (V.Ivantchenko)
+// 08-04-05 Major optimisation of internal interfaces (V.Ivanchenko)
+// 03-08-05 Angular correlations according to PRM (V.Ivanchenko)
 // 13-02-06 add ComputeCrossSectionPerAtom (mma)
 // 21-03-06 Fix problem of initialisation in case when cuts are not defined (VI)
 // 07-11-07 Improve sampling of final state (A.Bogdanov)
+// 28-02-08 Use precomputed Z^1/3 and Log(A) (V.Ivanchenko)
 //
 
@@ -74,12 +75,4 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-// static members
-//
-G4double G4MuBremsstrahlungModel::zdat[]={1., 4., 13., 29., 92.};
-G4double G4MuBremsstrahlungModel::adat[]={1.01, 9.01, 26.98, 63.55, 238.03};
-G4double G4MuBremsstrahlungModel::tdat[]={1.e3, 1.e4, 1.e5, 1.e6, 1.e7, 1.e8,
-                                          1.e9, 1.e10};
-
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -90,14 +83,15 @@
   : G4VEmModel(nam),
     particle(0),
+    sqrte(sqrt(exp(1.))),
+    bh(202.4),
+    bh1(446.),
+    btf(183.),
+    btf1(1429.),
+    fParticleChange(0),
     lowestKinEnergy(1.0*GeV),
-    minThreshold(1.0*keV),
-    nzdat(5),
-    ntdat(8),
-    NBIN(1000),
-    cutFixed(0.98*keV),
-    ignoreCut(false),
-    samplingTablesAreFilled(false)
+    minThreshold(1.0*keV)
 {
   theGamma = G4Gamma::Gamma();
+  nist = G4NistManager::Instance();
   if(p) SetParticle(p);
 }
@@ -125,14 +119,4 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4MuBremsstrahlungModel::SetParticle(const G4ParticleDefinition* p)
-{
-  if(!particle) {
-    particle = p;
-    mass = particle->GetPDGMass();
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
 void G4MuBremsstrahlungModel::Initialise(const G4ParticleDefinition* p,
                                          const G4DataVector& cuts)
@@ -142,4 +126,5 @@
   highKinEnergy = HighEnergyLimit();
 
+  // partial cross section is computed for fixed energy
   G4double fixedEnergy = 0.5*highKinEnergy;
 
@@ -148,19 +133,23 @@
   if(theCoupleTable) {
     G4int numOfCouples = theCoupleTable->GetTableSize();
-    
+
+    // clear old data    
     G4int nn = partialSumSigma.size();
     G4int nc = cuts.size();
     if(nn > 0) {
       for (G4int ii=0; ii<nn; ii++){
-	G4DataVector* a=partialSumSigma[ii];
+	G4DataVector* a = partialSumSigma[ii];
 	if ( a )  delete a;    
       } 
       partialSumSigma.clear();
     }
+    // fill new data
     if (numOfCouples>0) {
       for (G4int i=0; i<numOfCouples; i++) {
         G4double cute = DBL_MAX;
+
+	// protection for usage with extrapolator
         if(i < nc) cute = cuts[i];
-        if(cute < cutFixed || ignoreCut) cute = cutFixed;
+
         const G4MaterialCutsCouple* couple = 
 	  theCoupleTable->GetMaterialCutsCouple(i);
@@ -171,10 +160,13 @@
     }
   }
-  if(!samplingTablesAreFilled) MakeSamplingTables();
-  if(pParticleChange)
-    fParticleChange = 
-      reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-  else
-    fParticleChange = new G4ParticleChangeForLoss();
+
+  // define pointer to G4ParticleChange
+  if(!fParticleChange) {
+    if(pParticleChange)
+      fParticleChange = 
+	reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
+    else
+      fParticleChange = new G4ParticleChangeForLoss();
+  }
 }
 
@@ -188,24 +180,23 @@
 {
   G4double dedx = 0.0;
-  if (kineticEnergy <= lowestKinEnergy || ignoreCut) return dedx;
+  if (kineticEnergy <= lowestKinEnergy) return dedx;
 
   G4double tmax = kineticEnergy;
-  G4double cut  = min(cutEnergy,tmax);
-  if(cut < cutFixed) cut = cutFixed;
+  G4double cut  = std::min(cutEnergy,tmax);
+  if(cut < minThreshold) cut = minThreshold;
 
   const G4ElementVector* theElementVector = material->GetElementVector();
   const G4double* theAtomicNumDensityVector =
-                                          material->GetAtomicNumDensityVector();
+    material->GetAtomicNumDensityVector();
 
   //  loop for elements in the material
   for (size_t i=0; i<material->GetNumberOfElements(); i++) {
 
-    G4double Z = (*theElementVector)[i]->GetZ();
-    G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
-
-    G4double loss = ComputMuBremLoss(Z, A, kineticEnergy, cut);
+    G4double loss = 
+      ComputMuBremLoss((*theElementVector)[i]->GetZ(), kineticEnergy, cut);
 
     dedx += loss*theAtomicNumDensityVector[i];
   }
+  //  G4cout << "BR e= " << kineticEnergy << "  dedx= " << dedx << G4endl;
   if(dedx < 0.) dedx = 0.;
   return dedx;
@@ -214,5 +205,5 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z, G4double A,
+G4double G4MuBremsstrahlungModel::ComputMuBremLoss(G4double Z,
                                                    G4double tkin, G4double cut)
 {
@@ -239,5 +230,5 @@
     {
       G4double ep = (aa + xgi[i]*hhh)*totalEnergy;
-      loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
+      loss += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
     }
     aa += hhh;
@@ -254,5 +245,4 @@
                                            G4double tkin,
                                            G4double Z,
-                                           G4double A,
                                            G4double cut)
 {
@@ -272,5 +262,5 @@
   G4double bbb = log(vmax);
   G4int    kkk = (G4int)((bbb-aaa)/ak1)+k2 ;
-  G4double hhh = (bbb-aaa)/float(kkk);
+  G4double hhh = (bbb-aaa)/G4double(kkk);
 
   G4double aa = aaa;
@@ -281,5 +271,5 @@
     {
       G4double ep = exp(aa + xgi[i]*hhh)*totalEnergy;
-      cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, A, ep);
+      cross += ep*wgi[i]*ComputeDMicroscopicCrossSection(tkin, Z, ep);
     }
     aa += hhh;
@@ -287,4 +277,6 @@
 
   cross *=hhh;
+
+  //G4cout << "BR e= " << tkin<< "  cross= " << cross/barn << G4endl;
 
   return cross;
@@ -296,14 +288,7 @@
                                            G4double tkin,
                                            G4double Z,
-                                           G4double A,
                                            G4double gammaEnergy)
 //  differential cross section
 {
-  static const G4double sqrte=sqrt(exp(1.)) ;
-  static const G4double bh=202.4,bh1=446.,btf=183.,btf1=1429. ;
-  static const G4double rmass=mass/electron_mass_c2 ;
-  static const G4double cc=classic_electr_radius/rmass ;
-  static const G4double coeff= 16.*fine_structure_const*cc*cc/3. ;
-
   G4double dxsection = 0.;
 
@@ -315,20 +300,23 @@
   G4double rab0=delta*sqrte ;
 
-  G4double z13 = exp(-log(Z)/3.) ;
-  G4double dn  = 1.54*exp(0.27*log(A)) ;
+  G4int iz = G4int(Z);
+  if(iz < 1) iz = 1;
+
+  G4double z13 = 1.0/nist->GetZ13(iz);
+  G4double dn  = 1.54*nist->GetA27(iz);
 
   G4double b,b1,dnstar ;
 
-  if(Z<1.5)
+  if(1 == iz)
   {
-    b=bh;
-    b1=bh1;
-    dnstar=dn ;
+    b  = bh;
+    b1 = bh1;
+    dnstar = dn;
   }
   else
   {
-    b=btf;
-    b1=btf1;
-    dnstar = exp((1.-1./Z)*log(dn)) ;
+    b  = btf;
+    b1 = btf1;
+    dnstar = dn/std::pow(dn, 1./Z);
   }
 
@@ -359,48 +347,19 @@
                                            const G4ParticleDefinition*,
                                                  G4double kineticEnergy,
-						 G4double Z, G4double A,
+						 G4double Z, G4double,
                                                  G4double cutEnergy,
-                                                 G4double)
-{
-  G4double cut  = min(cutEnergy, kineticEnergy);
-  if(cut < cutFixed || ignoreCut) cut = cutFixed;
-  G4double cross =
-    ComputeMicroscopicCrossSection (kineticEnergy, Z, A/(g/mole), cut);
-  return cross;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuBremsstrahlungModel::CrossSectionPerVolume(
-					       const G4Material* material,
-                                               const G4ParticleDefinition*,
-                                                     G4double kineticEnergy,
-                                                     G4double cutEnergy,
-                                                     G4double maxEnergy)
+                                                 G4double maxEnergy)
 {
   G4double cross = 0.0;
-  if (cutEnergy >= maxEnergy || kineticEnergy <= lowestKinEnergy) return cross;
-  
-  G4double tmax = min(maxEnergy, kineticEnergy);
-  G4double cut  = min(cutEnergy, tmax);
-  if(cut < cutFixed || ignoreCut) cut = cutFixed;
-
-  const G4ElementVector* theElementVector = material->GetElementVector();
-  const G4double* theAtomNumDensityVector =
-                                          material->GetAtomicNumDensityVector();
-
-  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
-
-    G4double Z = (*theElementVector)[i]->GetZ();
-    G4double A = (*theElementVector)[i]->GetA()/(g/mole);
-
-    G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
-
-    if(tmax < kineticEnergy) {
-      cr -= ComputeMicroscopicCrossSection(kineticEnergy, Z, A, tmax);
-    }
-    cross += theAtomNumDensityVector[i] * cr;
-  }
-
+  if (kineticEnergy <= lowestKinEnergy) return cross;
+  G4double tmax = std::min(maxEnergy, kineticEnergy);
+  G4double cut  = std::min(cutEnergy, kineticEnergy);
+  if(cut < minThreshold) cut = minThreshold;
+  if (cut >= tmax) return cross;
+
+  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
+  if(tmax < kineticEnergy) {
+    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
+  }
   return cross;
 }
@@ -410,14 +369,15 @@
 G4DataVector* G4MuBremsstrahlungModel::ComputePartialSumSigma(
                                        const G4Material* material,
-                                             G4double kineticEnergy,
-                                             G4double cut)
-
-// Build the table of cross section per element. The table is built for MATERIAL
-// This table is used by DoIt to select randomly an element in the material.
+				       G4double kineticEnergy,
+				       G4double cut)
+
+// Build the table of cross section per element. 
+// The table is built for material 
+// This table is used to select randomly an element in the material.
 {
   G4int nElements = material->GetNumberOfElements();
   const G4ElementVector* theElementVector = material->GetElementVector();
   const G4double* theAtomNumDensityVector = 
-                                          material->GetAtomicNumDensityVector();
+    material->GetAtomicNumDensityVector();
 
   G4DataVector* dv = new G4DataVector();
@@ -426,9 +386,7 @@
 
   for (G4int i=0; i<nElements; i++ ) {
-
-    G4double Z = (*theElementVector)[i]->GetZ();
-    G4double A = (*theElementVector)[i]->GetA()/(g/mole) ;
     cross += theAtomNumDensityVector[i] 
-             * ComputeMicroscopicCrossSection(kineticEnergy, Z, A, cut);
+      * ComputeMicroscopicCrossSection(kineticEnergy, 
+				       (*theElementVector)[i]->GetZ(), cut);
     dv->push_back(cross);
   }
@@ -438,83 +396,18 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4MuBremsstrahlungModel::MakeSamplingTables()
-{
-
-  G4double AtomicNumber,AtomicWeight,KineticEnergy,
-           TotalEnergy,Maxep;
-
-  for (G4int iz=0; iz<nzdat; iz++)
-   {
-     AtomicNumber = zdat[iz];
-     AtomicWeight = adat[iz]*g/mole ;
-
-     for (G4int it=0; it<ntdat; it++)
-     {
-       KineticEnergy = tdat[it];
-       TotalEnergy = KineticEnergy + mass;
-       Maxep = KineticEnergy ;
-
-       G4double CrossSection = 0.0 ;
-
-       // calculate the differential cross section
-       // numerical integration in
-       //  log ...............
-       G4double c = log(Maxep/cutFixed) ;
-       G4double ymin = -5. ;
-       G4double ymax = 0. ;
-       G4double dy = (ymax-ymin)/NBIN ;
-
-       G4double y = ymin - 0.5*dy ;
-       G4double yy = ymin - dy ;
-       G4double x = exp(y);
-       G4double fac = exp(dy);
-       G4double dx = exp(yy)*(fac - 1.0);
-
-       for (G4int i=0 ; i<NBIN; i++)
-       {
-         y += dy ;
-         x *= fac;
-         dx*= fac;
-         G4double ep = cutFixed*exp(c*x) ;
-
-         CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
-                                                 KineticEnergy,AtomicNumber,
-                                                 AtomicWeight,ep) ;
-         ya[i]=y ;
-         proba[iz][it][i] = CrossSection ;
-
-       }
-
-       proba[iz][it][NBIN] = CrossSection ;
-       ya[NBIN] = 0. ;   //   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-       if(CrossSection > 0.)
-       {
-         for(G4int ib=0; ib<=NBIN; ib++)
-         {
-           proba[iz][it][ib] /= CrossSection ;
-         }
-       }
-     }
-   }
-  samplingTablesAreFilled = true;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MuBremsstrahlungModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
-						const G4MaterialCutsCouple* couple,
-						const G4DynamicParticle* dp,
-						G4double minEnergy,
-						G4double maxEnergy)
+void G4MuBremsstrahlungModel::SampleSecondaries(
+                              std::vector<G4DynamicParticle*>* vdp,
+			      const G4MaterialCutsCouple* couple,
+			      const G4DynamicParticle* dp,
+			      G4double minEnergy,
+			      G4double maxEnergy)
 {
   G4double kineticEnergy = dp->GetKineticEnergy();
   // check against insufficient energy
-  G4double tmax = min(kineticEnergy, maxEnergy);
-  G4double tmin = min(kineticEnergy, minEnergy);
-  if(tmin < cutFixed || ignoreCut) tmin = cutFixed;
+  G4double tmax = std::min(kineticEnergy, maxEnergy);
+  G4double tmin = std::min(kineticEnergy, minEnergy);
+  if(tmin < minThreshold) tmin = minThreshold;
   if(tmin >= tmax) return;
 
-  // ===== the begining of a new code  ======
   // ===== sampling of energy transfer ======
 
@@ -523,39 +416,30 @@
   // select randomly one element constituing the material
   const G4Element* anElement = SelectRandomAtom(couple);
+  G4double Z = anElement->GetZ();
 
   G4double totalEnergy   = kineticEnergy + mass;
   G4double totalMomentum = sqrt(kineticEnergy*(kineticEnergy + 2.0*mass));
 
-  G4double AtomicNumber = anElement->GetZ();
-  G4double AtomicWeight = anElement->GetA()/(g/mole);
-
-  G4double func1 = tmin*ComputeDMicroscopicCrossSection(
-				    kineticEnergy,AtomicNumber,
-				    AtomicWeight,tmin);
+  G4double func1 = tmin*
+    ComputeDMicroscopicCrossSection(kineticEnergy,Z,tmin);
 
   G4double lnepksi, epksi;
   G4double func2;
-  G4double ksi2;
 
   do {
     lnepksi = log(tmin) + G4UniformRand()*log(kineticEnergy/tmin);
     epksi   = exp(lnepksi);
-    func2   = epksi*ComputeDMicroscopicCrossSection(
-				kineticEnergy,AtomicNumber,
-				AtomicWeight,epksi);
-    ksi2 = G4UniformRand();
-
-  } while(func2/func1 < ksi2);
-
-  // ===== the end of a new code =====
-
-  // create G4DynamicParticle object for the Gamma
+    func2   = epksi*ComputeDMicroscopicCrossSection(kineticEnergy,Z,epksi);
+
+  } while(func2 < func1*G4UniformRand());
+
   G4double gEnergy = epksi;
 
-  // sample angle
+  // ===== sample angle =====
+
   G4double gam  = totalEnergy/mass;
-  G4double rmax = gam*min(1.0, totalEnergy/gEnergy - 1.0);
-  rmax *= rmax;
-  G4double x = G4UniformRand()*rmax/(1.0 + rmax);
+  G4double rmax = gam*std::min(1.0, totalEnergy/gEnergy - 1.0);
+  G4double rmax2= rmax*rmax;
+  G4double x = G4UniformRand()*rmax2/(1.0 + rmax2);
 
   G4double theta = sqrt(x/(1.0 - x))/gam;
@@ -577,5 +461,6 @@
 
   // save secondary
-  G4DynamicParticle* aGamma = new G4DynamicParticle(theGamma,gDirection,gEnergy);
+  G4DynamicParticle* aGamma = 
+    new G4DynamicParticle(theGamma,gDirection,gEnergy);
   vdp->push_back(aGamma);
 }
Index: trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuIonisation.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuIonisation.cc,v 1.54 2007/05/22 17:35:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuIonisation.cc,v 1.59 2009/02/26 11:04:20 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -86,4 +86,5 @@
 #include "G4MuBetheBlochModel.hh"
 #include "G4UniversalFluctuation.hh"
+#include "G4IonFluctuations.hh"
 #include "G4BohrFluctuations.hh"
 #include "G4UnitsTable.hh"
@@ -99,7 +100,8 @@
     isInitialised(false)
 {
-  SetStepFunction(0.2, 1*mm);
-  SetIntegral(true);
-  SetVerboseLevel(1);
+  //  SetStepFunction(0.2, 1*mm);
+  //SetIntegral(true);
+  //SetVerboseLevel(1);
+  SetProcessSubType(fIonisation);
 }
 
@@ -108,4 +110,22 @@
 G4MuIonisation::~G4MuIonisation()
 {}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4bool G4MuIonisation::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4MuIonisation::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....
@@ -122,18 +142,24 @@
     SetSecondaryParticle(G4Electron::Electron());
 
-    flucModel = new G4UniversalFluctuation();
+    // Bragg peak model
+    if (!EmModel(1)) SetEmModel(new G4BraggModel(),1);
+    EmModel(1)->SetLowEnergyLimit(MinKinEnergy());
+    EmModel(1)->SetHighEnergyLimit(0.2*MeV);
+    AddEmModel(1, EmModel(1), new G4IonFluctuations());
 
-    G4VEmModel* em = new G4BraggModel();
-    em->SetLowEnergyLimit(0.1*keV);
-    em->SetHighEnergyLimit(0.2*MeV);
-    AddEmModel(1, em, flucModel);
-    G4VEmModel* em1 = new G4BetheBlochModel();
-    em1->SetLowEnergyLimit(0.2*MeV);
-    em1->SetHighEnergyLimit(1.0*GeV);
-    AddEmModel(2, em1, flucModel);
-    G4VEmModel* em2 = new G4MuBetheBlochModel();
-    em2->SetLowEnergyLimit(1.0*GeV);
-    em2->SetHighEnergyLimit(100.0*TeV);
-    AddEmModel(3, em2, flucModel);
+    // high energy fluctuation model
+    if (!FluctModel()) SetFluctModel(new G4UniversalFluctuation());
+
+    // moderate energy model
+    if (!EmModel(2)) SetEmModel(new G4BetheBlochModel(),2);
+    EmModel(2)->SetLowEnergyLimit(0.2*MeV);
+    EmModel(2)->SetHighEnergyLimit(1.0*GeV);
+    AddEmModel(2, EmModel(2), FluctModel());
+
+    // high energy model
+    if (!EmModel(3)) SetEmModel(new G4MuBetheBlochModel(),3);
+    EmModel(3)->SetLowEnergyLimit(1.0*GeV);
+    EmModel(3)->SetHighEnergyLimit(MaxKinEnergy());
+    AddEmModel(3, EmModel(3), FluctModel());
 
     ratio = electron_mass_c2/mass;
@@ -145,10 +171,5 @@
 
 void G4MuIonisation::PrintInfo()
-{
-  G4cout << "      Bether-Bloch model for E > 0.2 MeV, "
-         << "parametrisation of Bragg peak below, "
-         << G4endl;
-  G4cout << "      radiative corrections for E > 1 GeV" << G4endl;
-}
+{}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/source/processes/electromagnetic/muons/src/G4MuMultipleScattering.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuMultipleScattering.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuMultipleScattering.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuMultipleScattering.cc,v 1.3 2007/11/09 19:48:10 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuMultipleScattering.cc,v 1.12 2008/10/16 13:37:04 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -----------------------------------------------------------------------------
@@ -47,5 +47,5 @@
 
 #include "G4MuMultipleScattering.hh"
-#include "G4MuMscModel.hh"
+#include "G4WentzelVIModel.hh"
 #include "G4MscStepLimitType.hh"
 
@@ -61,5 +61,5 @@
   samplez           = false ; 
   isInitialized     = false;  
-  SetRangeFactor(0.04);
+  SetRangeFactor(0.2);
   SetLateralDisplasmentFlag(true);
 }
@@ -84,23 +84,27 @@
   if(isInitialized) {
 
-    if (p->GetParticleType() != "nucleus") {
+    if (p->GetParticleType() != "nucleus" && p->GetPDGMass() < GeV) {
       mscModel->SetStepLimitType(StepLimitType());
       mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
-      //mscModel->SetThetaLimit(thetaLimit);
       mscModel->SetRangeFactor(RangeFactor());
     }
+    mscModel->SetPolarAngleLimit(PolarAngleLimit());
     return;
   }
 
-  if (p->GetParticleType() == "nucleus") {
+  if (p->GetParticleType() == "nucleus" || p->GetPDGMass() > GeV) {
     SetLateralDisplasmentFlag(false);
     SetBuildLambdaTable(false);
-    //    SetRangeFactor(0.2);
   }
 
-  // initialisation of parameters
-  //  G4String part_name = p->GetParticleName();
-  mscModel = new G4MuMscModel(RangeFactor(),thetaLimit);
+  // initialisation of the model
+
+  mscModel = new G4WentzelVIModel();
+  mscModel->SetStepLimitType(StepLimitType());
   mscModel->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
+  mscModel->SetRangeFactor(RangeFactor());
+  mscModel->SetPolarAngleLimit(PolarAngleLimit());
+  mscModel->SetLowEnergyLimit(MinKinEnergy());
+  mscModel->SetHighEnergyLimit(MaxKinEnergy());
 
   AddEmModel(1,mscModel);
@@ -112,8 +116,8 @@
 void G4MuMultipleScattering::PrintInfo()
 {
-  G4cout << "      Boundary/stepping algorithm is active with RangeFactor= "
-	 << RangeFactor()
-	 << "  Step limit type " << StepLimitType()
-	 << G4endl;
+  G4cout << "      RangeFactor= " << RangeFactor()
+         << ", step limit type: " << StepLimitType()
+         << ", lateralDisplacement: " << LateralDisplasmentFlag()
+         << G4endl;
 }
 
Index: trunk/source/processes/electromagnetic/muons/src/G4MuPairProduction.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuPairProduction.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuPairProduction.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuPairProduction.cc,v 1.48 2007/05/22 17:35:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuPairProduction.cc,v 1.52 2009/02/20 14:48:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -89,5 +89,7 @@
     lowestKinEnergy(1.*GeV),
     isInitialised(false)
-{}
+{
+  SetProcessSubType(fPairProdByCharged);
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -98,6 +100,23 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4MuPairProduction::InitialiseEnergyLossProcess(const G4ParticleDefinition* part,
-                                                     const G4ParticleDefinition*)
+G4bool G4MuPairProduction::IsApplicable(const G4ParticleDefinition& p)
+{
+  return (p.GetPDGCharge() != 0.0 && p.GetPDGMass() > 10.0*MeV);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4MuPairProduction::MinPrimaryEnergy(const G4ParticleDefinition*,
+					      const G4Material*,
+					      G4double)
+{
+  return lowestKinEnergy;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4MuPairProduction::InitialiseEnergyLossProcess(
+                         const G4ParticleDefinition* part,
+			 const G4ParticleDefinition*)
 {
   if (!isInitialised) {
@@ -110,7 +129,7 @@
     G4MuPairProductionModel* em = new G4MuPairProductionModel();
     em->SetLowestKineticEnergy(lowestKinEnergy);
-    G4VEmFluctuationModel* fm = new G4UniversalFluctuation();
-    em->SetLowEnergyLimit(0.1*keV);
-    em->SetHighEnergyLimit(100.0*TeV);
+    G4VEmFluctuationModel* fm = 0;
+    em->SetLowEnergyLimit(MinKinEnergy());
+    em->SetHighEnergyLimit(MaxKinEnergy());
     AddEmModel(1, em, fm);
   }
@@ -120,8 +139,5 @@
 
 void G4MuPairProduction::PrintInfo()
-{
-  G4cout << "      Parametrised model "
-         << G4endl;
-}
+{}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
Index: trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/muons/src/G4MuPairProductionModel.cc	(revision 961)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4MuPairProductionModel.cc,v 1.35 2007/10/11 13:52:04 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4MuPairProductionModel.cc,v 1.40 2009/02/20 14:48:16 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -104,11 +104,12 @@
                                                  const G4String& nam)
   : G4VEmModel(nam),
-  minPairEnergy(4.*electron_mass_c2),
-  lowestKinEnergy(1.*GeV),
-  factorForCross(4.*fine_structure_const*fine_structure_const
+    particle(0),
+    factorForCross(4.*fine_structure_const*fine_structure_const
                    *classic_electr_radius*classic_electr_radius/(3.*pi)),
     sqrte(sqrt(exp(1.))),
     currentZ(0),
-    particle(0),
+    fParticleChange(0),
+    minPairEnergy(4.*electron_mass_c2),
+    lowestKinEnergy(1.*GeV),
     nzdat(5),
     ntdat(8),
@@ -118,8 +119,8 @@
     ymax(0.),
     dy((ymax-ymin)/nbiny),
-    ignoreCut(false),
     samplingTablesAreFilled(false)
 {
   SetLowEnergyLimit(minPairEnergy);
+  nist = G4NistManager::Instance();
 
   theElectron = G4Electron::Electron();
@@ -144,10 +145,9 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4MuPairProductionModel::SetParticle(const G4ParticleDefinition* p)
-{
-  if(!particle) {
-    particle = p;
-    particleMass = particle->GetPDGMass();
-  }
+G4double G4MuPairProductionModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
+						     G4double kineticEnergy)
+{
+  G4double maxPairEnergy = kineticEnergy + particleMass*(1.0 - 0.75*sqrte*z13);
+  return maxPairEnergy;
 }
 
@@ -161,17 +161,10 @@
     MakeSamplingTables();
   }
-  if(pParticleChange) {
-    if(ignoreCut) {
-      gParticleChange = 
-	reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
-      fParticleChange = 0;
-    } else {
+  if(!fParticleChange) {
+    if(pParticleChange) 
       fParticleChange = 
 	reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
-      gParticleChange = 0;
-    }
-  } else {
-    fParticleChange = new G4ParticleChangeForLoss();
-    gParticleChange = 0;
+    else 
+      fParticleChange = new G4ParticleChangeForLoss();
   }
 }
@@ -186,6 +179,5 @@
 {
   G4double dedx = 0.0;
-  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy 
-      || ignoreCut)
+  if (cutEnergy <= minPairEnergy || kineticEnergy <= lowestKinEnergy)
     return dedx;
 
@@ -216,5 +208,5 @@
   G4double loss = 0.0;
 
-  G4double cut  = min(cutEnergy,tmax);
+  G4double cut = std::min(cutEnergy,tmax);
   if(cut <= minPairEnergy) return loss;
 
@@ -251,11 +243,8 @@
                                            G4double Z,
                                            G4double cut)
-
-{
-  G4double cross = 0. ;
-
+{
+  G4double cross = 0.;
   SetCurrentElement(Z);
   G4double tmax = MaxSecondaryEnergy(particle, tkin);
-
   if (tmax <= cut) return cross;
 
@@ -400,42 +389,18 @@
 						 G4double Z, G4double,
                                                  G4double cutEnergy,
-                                                 G4double)
-{
-  G4double cut  = max(minPairEnergy,cutEnergy);
-  if(ignoreCut) cut = minPairEnergy;
-  G4double cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
-  return cross;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuPairProductionModel::CrossSectionPerVolume(
-					       const G4Material* material,
-                                               const G4ParticleDefinition*,
-                                                     G4double kineticEnergy,
-                                                     G4double cutEnergy,
-						     G4double maxEnergy)
+                                                 G4double maxEnergy)
 {
   G4double cross = 0.0;
   if (kineticEnergy <= lowestKinEnergy) return cross;
 
-  maxEnergy += particleMass;
-
-  const G4ElementVector* theElementVector = material->GetElementVector();
-  const G4double* theAtomNumDensityVector = material->
-                                                    GetAtomicNumDensityVector();
-
-  for (size_t i=0; i<material->GetNumberOfElements(); i++) {
-    G4double Z = (*theElementVector)[i]->GetZ();
-    SetCurrentElement(Z);
-    G4double tmax = min(maxEnergy,MaxSecondaryEnergy(particle, kineticEnergy));
-    G4double cut  = max(minPairEnergy,cutEnergy);
-    if(ignoreCut) cut = minPairEnergy;
-    if(cut < tmax) {
-      G4double cr = ComputeMicroscopicCrossSection(kineticEnergy, Z, cut)
-                  - ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
-
-      cross += theAtomNumDensityVector[i] * cr;
-    }
+  SetCurrentElement(Z);
+  G4double tmax = std::min(maxEnergy, kineticEnergy);
+  G4double cut  = std::min(cutEnergy, kineticEnergy);
+  if(cut < minPairEnergy) cut = minPairEnergy;
+  if (cut >= tmax) return cross;
+
+  cross = ComputeMicroscopicCrossSection (kineticEnergy, Z, cut);
+  if(tmax < kineticEnergy) {
+    cross -= ComputeMicroscopicCrossSection(kineticEnergy, Z, tmax);
   }
   return cross;
@@ -451,35 +416,42 @@
     SetCurrentElement(Z);
 
-    for (G4int it=0; it<ntdat; it++)
-    {
+    for (G4int it=0; it<ntdat; it++) {
+
       G4double kineticEnergy = tdat[it];
       G4double maxPairEnergy = MaxSecondaryEnergy(particle,kineticEnergy);
-
+      // G4cout << "Z= " << currentZ << " z13= " << z13 
+      //<< " mE= " << maxPairEnergy << G4endl;
       G4double CrossSection = 0.0 ;
 
-      G4double y = ymin - 0.5*dy ;
-      G4double yy = ymin - dy ;
-      G4double x = exp(y);
-      G4double fac = exp(dy);
-      G4double dx = exp(yy)*(fac - 1.0);
-
-      G4double c = log(maxPairEnergy/minPairEnergy);
-
-      for (G4int i=0 ; i<nbiny; i++)
-      {
-        y += dy ;
-        if(c > 0.0) {
-	  x *= fac;
-          dx*= fac;
-          G4double ep = minPairEnergy*exp(c*x) ;
-          CrossSection += ep*dx*ComputeDMicroscopicCrossSection(
-                                      kineticEnergy, Z, ep);
-        }
-        ya[i] = y;
-        proba[iz][it][i] = CrossSection;
+      if(maxPairEnergy > minPairEnergy) {
+
+	G4double y = ymin - 0.5*dy ;
+	G4double yy = ymin - dy ;
+	G4double x = exp(y);
+	G4double fac = exp(dy);
+	G4double dx = exp(yy)*(fac - 1.0);
+
+	G4double c = log(maxPairEnergy/minPairEnergy);
+
+	for (G4int i=0 ; i<nbiny; i++) {
+	  y += dy ;
+	  if(c > 0.0) {
+	    x *= fac;
+	    dx*= fac;
+	    G4double ep = minPairEnergy*exp(c*x) ;
+	    CrossSection += 
+	      ep*dx*ComputeDMicroscopicCrossSection(kineticEnergy, Z, ep);
+	  }
+	  ya[i] = y;
+	  proba[iz][it][i] = CrossSection;
+	}
+       
+      } else {
+	for (G4int i=0 ; i<nbiny; i++) {
+	  proba[iz][it][i] = CrossSection;
+	}
       }
 
       ya[nbiny]=ymax;
-
       proba[iz][it][nbiny] = CrossSection;
 
@@ -515,5 +487,5 @@
   G4double maxEnergy     = std::min(tmax, maxPairEnergy);
   G4double minEnergy     = std::max(tmin, minPairEnergy);
-  if(ignoreCut)minEnergy = minPairEnergy;
+
   if(minEnergy >= maxEnergy) return;
   //G4cout << "emin= " << minEnergy << " emax= " << maxEnergy 
@@ -618,8 +590,5 @@
   // primary change
   kineticEnergy -= (ElectronEnergy + PositronEnergy);
-  if(fParticleChange)
-    fParticleChange->SetProposedKineticEnergy(kineticEnergy);
-  else 
-    gParticleChange->SetProposedKineticEnergy(kineticEnergy);
+  fParticleChange->SetProposedKineticEnergy(kineticEnergy);
 
   vdp->push_back(aParticle1);
@@ -655,5 +624,4 @@
     G4double maxPairEnergy = MaxSecondaryEnergy(particle,kinEnergy);
     G4double minEnergy     = std::max(tmin, minPairEnergy);
-    if(ignoreCut)minEnergy = minPairEnergy;
 
     G4int iz;
