Index: trunk/source/processes/electromagnetic/muons/include/G4MuMscModel.hh
===================================================================
--- trunk/source/processes/electromagnetic/muons/include/G4MuMscModel.hh	(revision 1007)
+++ 	(revision )
@@ -1,256 +1,0 @@
-//
-// ********************************************************************
-// * 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: G4MuMscModel.hh,v 1.4 2007/11/09 19:48:09 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
-//
-// -------------------------------------------------------------------
-//
-//
-// GEANT4 Class header file
-//
-//
-// File name:     G4MuMscModel
-//
-// Author:        V.Ivanchenko on base of L.Urban model
-//
-// Creation date: 25.10.2007
-//
-// Modifications:
-//
-//
-// Class Description:
-//
-// Implementation of the model of multiple scattering based on
-// H.W.Lewis Phys Rev 78 (1950) 526 and L.Urban model
-
-// -------------------------------------------------------------------
-//
-
-#ifndef G4MuMscModel_h
-#define G4MuMscModel_h 1
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-#include "G4eCoulombScatteringModel.hh"
-#include "G4PhysicsTable.hh"
-#include "G4MscStepLimitType.hh"
-#include "G4MaterialCutsCouple.hh"
-
-class G4LossTableManager;
-class G4ParticleChangeForMSC;
-class G4SafetyHelper;
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-class G4MuMscModel : public G4eCoulombScatteringModel
-{
-
-public:
-
-  G4MuMscModel(G4double frange = 0.2,
-	       G4double thetaMax = 0.04,
-	       G4double tMax = TeV*TeV, 
-	       const G4String& nam = "MuMscUni");
-
-  virtual ~G4MuMscModel();
-
-  void Initialise(const G4ParticleDefinition*, const G4DataVector&);
-
-  G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
-				      G4double KineticEnergy,
-				      G4double AtomicNumber,
-				      G4double AtomicWeight=0., 
-				      G4double cut = DBL_MAX,
-				      G4double emax= DBL_MAX);
-
-  void SampleScattering(const G4DynamicParticle*, G4double safety);
-
-  void SampleSecondaries(std::vector<G4DynamicParticle*>*, 
-			 const G4MaterialCutsCouple*,
-			 const G4DynamicParticle*,
-			 G4double,
-			 G4double);
-
-  G4double ComputeTruePathLengthLimit(const G4Track& track,
-				      G4PhysicsTable* theLambdaTable,
-				      G4double currentMinimalStep);
-
-  G4double ComputeGeomPathLength(G4double truePathLength);
-
-  G4double ComputeTrueStepLength(G4double geomStepLength);
-
-  inline void SetStepLimitType(G4MscStepLimitType);
-
-  inline void SetLateralDisplasmentFlag(G4bool val);
-
-  inline G4double GetLambda(G4double kinEnergy);
-
-  inline G4double GetLambda2(G4double kinEnergy);
-
-  //  inline void SetThetaLimit(G4double);
-
-  inline void SetRangeFactor(G4double);
-
-private:
-
-  void BuildTables();
-
-  G4double ComputeLambda2(G4double kinEnergy, G4double cut);
-
-  inline void DefineMaterial(const G4MaterialCutsCouple*);
-
-  //  hide assignment operator
-  G4MuMscModel & operator=(const  G4MuMscModel &right);
-  G4MuMscModel(const  G4MuMscModel&);
-
-  G4ParticleChangeForMSC*     fParticleChange;
-
-  G4SafetyHelper*             safetyHelper;
-  G4PhysicsTable*             theLambdaTable;
-  G4PhysicsTable*             theLambda2Table;
-  G4LossTableManager*         theManager;
-  const G4DataVector*         currentCuts;
-
-  G4double dtrl;
-  G4double facrange;
-  G4double thetaLimit;
-  G4double numlimit;
-  G4double tlimitminfix;
-  G4double invsqrt12;
-  G4double lowBinEnergy;
-  G4double highBinEnergy;
-
-  // cash
-  G4double preKinEnergy;
-  G4double xSection;
-  G4double ecut;
-  G4double lambda0;
-  G4double tPathLength;
-  G4double zPathLength;
-  G4double lambdaeff;
-  G4double currentRange; 
-  G4double par1;
-  G4double par2;
-  G4double par3;
-
-  G4int    currentMaterialIndex;
-
-  G4int    nbins;
-  G4int    nwarnings;
-  G4int    nwarnlimit;
-
-  const G4MaterialCutsCouple* currentCouple;
-
-  G4MscStepLimitType steppingAlgorithm;
-
-  G4bool   samplez;
-  G4bool   latDisplasment;
-  G4bool   isInitialized;
-  G4bool   buildTables;
-  G4bool   newrun;
-  G4bool   inside;
-};
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-void G4MuMscModel::SetLateralDisplasmentFlag(G4bool val) 
-{ 
-  latDisplasment = val;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-/*
-inline
-void G4MuMscModel::SetThetaLimit(G4double val) 
-{ 
-  thetaLimit = val;
-}
-*/
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-void G4MuMscModel::SetRangeFactor(G4double val) 
-{ 
-  facrange = val;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-void G4MuMscModel::SetStepLimitType(G4MscStepLimitType val) 
-{ 
-  steppingAlgorithm = val;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-void G4MuMscModel::DefineMaterial(const G4MaterialCutsCouple* cup) 
-{ 
-  if(cup != currentCouple) {
-    currentCouple = cup;
-    currentMaterialIndex = currentCouple->GetIndex(); 
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-G4double G4MuMscModel::GetLambda(G4double e)
-{
-  G4double x;
-  if(theLambdaTable) {
-    G4bool b;
-    x = ((*theLambdaTable)[currentMaterialIndex])->GetValue(e, b);
-  } else {
-    x = CrossSection(currentCouple,particle,e);
-  }
-  if(x > DBL_MIN) x = 1./x;
-  else            x = DBL_MAX;
-  return x;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-inline
-G4double G4MuMscModel::GetLambda2(G4double e)
-{
-  G4double x;
-  if(theLambda2Table) {
-    G4bool b;
-    x = ((*theLambda2Table)[currentMaterialIndex])->GetValue(e, b);
-  } else {
-    x = ComputeLambda2(e, (*currentCuts)[currentMaterialIndex]);
-  }
-  return x;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-#endif
-
Index: trunk/source/processes/electromagnetic/muons/src/G4MuMscModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/muons/src/G4MuMscModel.cc	(revision 1007)
+++ 	(revision )
@@ -1,571 +1,0 @@
-//
-// ********************************************************************
-// * 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: G4MuMscModel.cc,v 1.6 2007/11/11 17:40:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
-//
-// -------------------------------------------------------------------
-//
-// GEANT4 Class file
-//
-//
-// File name:   G4MuMscModel
-//
-// Author:      Laszlo Mu
-//
-// Creation date: 03.03.2001
-//
-// Modifications:
-//
-// 27-03-03 Move model part from G4MultipleScattering80 (V.Ivanchenko)
-//
-
-// Class Description:
-//
-// Implementation of the model of multiple scattering based on
-// H.W.Lewis Phys Rev 78 (1950) 526 and others
-
-// -------------------------------------------------------------------
-//
-
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-#include "G4MuMscModel.hh"
-#include "Randomize.hh"
-#include "G4Electron.hh"
-#include "G4LossTableManager.hh"
-#include "G4ParticleChangeForMSC.hh"
-#include "G4TransportationManager.hh"
-#include "G4SafetyHelper.hh"
-#include "G4eCoulombScatteringModel.hh"
-#include "G4PhysicsTableHelper.hh"
-#include "G4ElementVector.hh"
-#include "G4ProductionCutsTable.hh"
-#include "G4PhysicsLogVector.hh"
-//#include "G4Poisson.hh"
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-using namespace std;
-
-G4MuMscModel::G4MuMscModel(G4double frange, 
-			   G4double thetaMax, 
-			   G4double tMax,  
-			   const G4String& nam)
-  : G4eCoulombScatteringModel(0.0,thetaMax,false,tMax,nam),
-    theLambdaTable(0),
-    theLambda2Table(0),
-    dtrl(0.05),
-    facrange(frange),
-    thetaLimit(thetaMax),
-    numlimit(0.2),
-    lowBinEnergy(keV),
-    highBinEnergy(PeV),
-    nbins(60),
-    nwarnings(0),
-    nwarnlimit(50),
-    currentCouple(0),
-    isInitialized(false),
-    buildTables(true),
-    newrun(true),
-    inside(false)
-{
-  invsqrt12 = 1./sqrt(12.);
-  tlimitminfix = 1.e-6*mm;
-  theManager = G4LossTableManager::Instance(); 
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4MuMscModel::~G4MuMscModel()
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MuMscModel::Initialise(const G4ParticleDefinition* p,
-			      const G4DataVector& cuts)
-{
-  SetupParticle(p);
-  newrun = true;
-  xSection = currentRange = targetZ = ecut = tkin = 0.0;
-  // set values of some data members
-  if(!isInitialized) {
-    isInitialized = true;
-    if(p->GetParticleName() == "GenericIon") buildTables = false;
-
-    if (pParticleChange)
-      fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
-    else
-      fParticleChange = new G4ParticleChangeForMSC();
-
-    safetyHelper = G4TransportationManager::GetTransportationManager()
-      ->GetSafetyHelper();
-    safetyHelper->InitialiseHelper();
-  }
-  G4eCoulombScatteringModel::Initialise(p, cuts);
-  currentCuts = &cuts;
-  if(buildTables)
-    theLambda2Table = G4PhysicsTableHelper::PreparePhysicsTable(theLambda2Table);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MuMscModel::BuildTables()
-{
-  //G4cout << "G4MuMscModel::BuildTables flags newrun= " << newrun 
-  //	 << "  buildTables= " << buildTables << G4endl;
-  newrun = false;
-  if(!buildTables) return;
-
-  // Access to materials
-  const G4ProductionCutsTable* theCoupleTable=
-        G4ProductionCutsTable::GetProductionCutsTable();
-  size_t numOfCouples = theCoupleTable->GetTableSize();
-  G4double e, s, cut;
-
-  for(size_t i=0; i<numOfCouples; i++) {
-
-    if (theLambda2Table->GetFlag(i)) {
-
-      // create physics vector and fill it
-      DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
-      cut = (*currentCuts)[currentMaterialIndex];
-      G4PhysicsVector* aVector =
-	new G4PhysicsLogVector(lowBinEnergy, highBinEnergy, nbins);
-      for(G4int j=0; j<nbins; j++) {
-        e = aVector->GetLowEdgeEnergy(j);
-        s = ComputeLambda2(e, cut);
-        //G4cout << j << "  " << currentCouple->GetMaterial()->GetName() 
-	//       << "  e(MeV)= " << e << " cut(MeV)= " << cut
-	//       << " L2= " << s << G4endl;  
-	aVector->PutValue(j, s);
-      }
-      
-      G4PhysicsTableHelper::SetPhysicsVector(theLambda2Table, i, aVector);
-    }
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuMscModel::ComputeCrossSectionPerAtom( 
-                             const G4ParticleDefinition* p,
-			     G4double kinEnergy,
-			     G4double Z, G4double A,
-			     G4double cutEnergy, G4double)
-{
-  if(p == particle && kinEnergy == tkin && Z == targetZ &&
-     cutEnergy == ecut) return xSection;
-  ecut = cutEnergy;
-  xSection = 0.0;
-  SetupParticle(p);
-  G4double ekin = std::max(keV, kinEnergy);
-  SetupTarget(Z, A, ekin);
-
-  G4double tmax = tkin;
-  if(p == theElectron) tmax *= 0.5;
-  else if(p != thePositron) {
-    G4double ratio = electron_mass_c2/mass;
-    tmax = 2.0*mom2/
-      (electron_mass_c2*(1.0 + ratio*(tkin/mass + 1.0) + ratio*ratio));
-  }
-  G4double t = std::min(cutEnergy, tmax);
-  G4double mom21 = t*(t + 2.0*electron_mass_c2);
-  t = tkin - t;
-  G4double mom22 = t*(t + 2.0*mass);
-  cosTetMaxElec = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
-  if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
-
-  if(cosTetMaxElec < 1.0) {
-    G4double x2 = screenZ/(1.0 - cosTetMaxElec + screenZ);
-    xSection += (x2 - 1.0 - log(x2))/Z;
-  }
-  //  G4cout << "cut= " << ecut << " e= " << tkin << " croosE= " 
-  //  << xSection/barn << G4endl;
-
-  if(cosTetMaxNuc < 1.0) {
-    G4double x1 = screenZ*formfactA;
-    G4double x2 = 1.0 - cosTetMaxNuc + screenZ;
-    G4double x3 = 1.0 - x1; 
-    G4double x4 = 1.0/(formfactA*x2 + x3);
-    G4double x5 = screenZ/x2;
-    xSection += ((1.0 - 2.0*x1/x3)*log(x4/x5) - 1.0 + 
-		 x5 - (1.0 - 4.0*x1)*(1.0 - x4))/(x3*x3); 
-  }
-  xSection *= coeff*Z*Z*chargeSquare*invbeta2/mom2; 
-  //  G4cout << " croosE= " << xSection/barn << " screenZ= " 
-  //	 << screenZ << " formF= " << formfactA << G4endl;
-  return xSection; 
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuMscModel::ComputeLambda2(G4double kinEnergy,
-				      G4double cutEnergy)
-{
-  G4double res = 0.0;
-  SetupParticle(particle);
-  G4double ekin = std::max(keV, kinEnergy);
-
-  const G4Material* mat = currentCouple->GetMaterial();
-  const G4ElementVector* theElementVector = mat->GetElementVector();
-  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
-  size_t nelm = mat->GetNumberOfElements();
-
-  SetupKinematic(ekin);
-
-  G4double tmax = tkin;
-  if(particle == theElectron) tmax *= 0.5;
-  else if(particle != thePositron) {
-    G4double ratio = electron_mass_c2/mass;
-    tmax = 2.0*mom2/
-      (electron_mass_c2*(1.0 + ratio*(tkin/mass + 1.0) + ratio*ratio));
-  }
-  G4double t = std::min(cutEnergy, tmax);
-  G4double mom21 = t*(t + 2.0*electron_mass_c2);
-  t = tkin - t;
-  G4double mom22 = t*(t + 2.0*mass);
-  cosTetMaxElec = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
-  if(cosTetMaxElec < 0.0) cosTetMaxElec = 0.0;
-
-  G4double x, x1, x2, y;
-
-  for (size_t i=0; i<nelm; i++) {
-    const G4Element* elm = (*theElementVector)[i];
-    G4double Z = elm->GetZ();
-    SetupTarget(Z, elm->GetN(), tkin);
-    G4double s = 0.0;  
-    G4double costm = cosTetMaxElec;
-    if(costm < cosTetMaxNuc) costm = cosTetMaxNuc;
-    if(costm < 1.0) {
-      x = 1.0 - costm + screenZ;
-      y = (x - screenZ*(screenZ/x + 2.0*log(x/screenZ)))/Z;
-      if(y < 0.0) {
-        nwarnings++;
-        if(nwarnings < nwarnlimit) 
-	  G4cout << "Electron scattering <0 for L2 " << y << G4endl;
-        y = 0.0;
-      }
-      s += y;
-    }
-    //  G4cout << "cut= " << cut << " e= " << tkin << " croosE= " 
-    //  << xSection/barn << G4endl;
-
-    // limit main integral because of nuclear size effect
-
-    if(cosTetMaxNuc < 1.0) {
-      x1 = screenZ*formfactA;
-      x2 = 1.0 - cosTetMaxNuc + screenZ;
-      G4double x3 = 1.0 - x1; 
-      G4double f  = 1.0/formfactA;
-      G4double d  = f - screenZ;
-      G4double x4 = f/(x2 + d);
-      G4double x5 = screenZ/x2;
-      y = (screenZ*(1.0 - x5) + (d*d - screenZ*(2.0*d - 3.0*screenZ))*(1.0 - x4)/f -
-	   2.0*screenZ*f*log(x4/x5)/d)/(x3*x3); 
-      if(y < 0.0) {
-        nwarnings++;
-        if(nwarnings < nwarnlimit) 
-	  G4cout << "Nuclear scattering <0 for L2 " << y << G4endl;
-        y = 0.0;
-      }
-      s += y;
-    }
-
-    res += Z*Z*s*theAtomNumDensityVector[i]; 
-  }
-  res *= 0.25*coeff*chargeSquare*invbeta2/mom2; 
-  //  G4cout << " croosE= " << xSection/barn << " screenZ= " 
-  //	 << screenZ << " formF= " << formfactA << G4endl;
-  return res; 
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuMscModel::ComputeTruePathLengthLimit(
-                             const G4Track& track,
-			     G4PhysicsTable* theTable,
-			     G4double currentMinimalStep)
-{
-  G4double tlimit = currentMinimalStep;
-  const G4DynamicParticle* dp = track.GetDynamicParticle();
-
-  // initialisation for 1st step  
-  if(track.GetCurrentStepNumber() == 1) {
-    inside = false;
-    SetupParticle(dp->GetDefinition());
-    theLambdaTable = theTable;
-    if(newrun && buildTables) BuildTables();
-  }
-
-  // initialisation for each step  
-  preKinEnergy = dp->GetKineticEnergy();
-  DefineMaterial(track.GetMaterialCutsCouple());
-  lambda0 = GetLambda(preKinEnergy);
-  currentRange = 
-    theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
-
-  // extra check for abnormal situation
-  // this check needed to run MSC with eIoni and eBrem inactivated
-  if(tlimit > currentRange) tlimit = currentRange;
-
-  // stop here if small range particle
-  if(inside) return tlimit;   
-
-  // pre step
-  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
-  G4StepStatus stepStatus = sp->GetStepStatus();
-  G4double presafety = sp->GetSafety();
-
-  // compute presafety again if presafety <= 0 and no boundary
-  // i.e. when it is needed for optimization purposes
-  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) 
-    presafety = safetyHelper->ComputeSafety(sp->GetPosition()); 
-
-  //  G4cout << "G4MuMscModel::ComputeTruePathLengthLimit tlimit= " 
-  //	 <<tlimit<<" safety= " << presafety
-  //	 << " range= " <<currentRange<<G4endl;
-
-  // far from geometry boundary
-  if(currentRange < presafety) {
-    inside = true;  
-    
-    // limit mean scattering angle
-  } else {
-    tlimit = std::min(facrange*lambda0, tlimit);
-  }
-  /*
-  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
-	 << " L0= " << lambda0 << " R= " << currentRange
-	 << "tlimit= " << tlimit  
-  	 << " currentMinimalStep= " << currentMinimalStep << G4endl;
-  */
-  return tlimit;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuMscModel::ComputeGeomPathLength(G4double truelength)
-{
-  tPathLength  = truelength;
-  zPathLength  = tPathLength;
-
-  G4double tau = tPathLength/lambda0;
-  lambdaeff    = lambda0;
-  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
-  //	 << " lambda0= " << lambda0 << " tau= " << tau << G4endl; 
-  // small step
-  if(tau < numlimit) {
-    par1 = -1. ;  
-    par2 = par3 = 0. ;  
-    zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
-
-    // medium step
-  } else if(tPathLength < currentRange*dtrl) {
-    zPathLength = lambda0*(1.0 - exp(-tau));
-
-  } else if(tkin < mass) {
-
-    par1 = 1./currentRange;
-    par2 = 1./(par1*lambda0);
-    par3 = 1.+ par2;
-    lambdaeff = 1.0/(par1*par3);
-    G4double x = tPathLength/currentRange;
-    G4double x1;
-    if(x < numlimit) x1 = x*(1.0  - 0.5*x + x*x/3.0);
-    else             x1 = log(1.0 - x); 
-    zPathLength = lambdaeff*(1.-exp(par3*x1));
-
-  } else {
-
-    G4double T1 = theManager->GetEnergy(particle,
-					currentRange-tPathLength,
-					currentCouple);
-    G4double lambda1 = GetLambda(T1);
-
-    par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
-    par2 = 1./(par1*lambda0) ;
-    par3 = 1.+ par2 ;
-    lambdaeff = 1.0/(par1*par3);
-    zPathLength = lambdaeff*(1.-exp(par3*log(lambda1/lambda0)));
-  }
-
-  //  if(zPathLength > lambda0) zPathLength = lambda0;
-  if(zPathLength > tPathLength) zPathLength = tPathLength;
-
-  return zPathLength;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4double G4MuMscModel::ComputeTrueStepLength(G4double geomStepLength)
-{
-  // step defined other than transportation 
-  if(geomStepLength == zPathLength) return tPathLength;
-
-  tPathLength  = geomStepLength;
-  zPathLength  = geomStepLength;
-  G4double tau = geomStepLength/lambda0;
-  if(tau < numlimit) {
-    tPathLength *= (1.0 + 0.5*tau - tau*tau/3.0); 
-  
-  } else if(par1 <  0.) {
-    tPathLength = -lambda0*log(1.0 - tau); 
-
-  } else {
-    G4double x = par1*par3*geomStepLength;
-    if(x < numlimit) 
-      tPathLength = (1.- exp(- x*(1.- 0.5*x + x*x/3.0)/par3))/par1 ;
-    else if (x < 1.0)  
-      tPathLength = (1.-exp(log(1.- x)/par3))/par1;
-    else 
-      tPathLength = currentRange;
-  }
-  if(tPathLength < geomStepLength) tPathLength = geomStepLength;
-
-  return tPathLength;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MuMscModel::SampleScattering(const G4DynamicParticle* dynParticle,
-				    G4double safety)
-{
-  G4double kinEnergy = dynParticle->GetKineticEnergy();
-  if(kinEnergy == 0.0) return;
-  G4double x1  = 0.5*tPathLength/lambdaeff;
-  
-  /*
-  G4cout << "G4MuMscModel::SampleScattering t(mm)= " << tPathLength
-	 << " 1/lambdaeff= " << 1.0/lambdaeff 
-	 << " matIdx= " << currentMaterialIndex << G4endl;
-  */
-  /*
-  G4double y1  = 1.0 - x1;
-  G4double x2  = tPathLength*GetLambda2(0.5*(preKinEnergy + kinEnergy)); 
-  G4double x3  = (x2 - x1*x1)/(x1*y1);
-  if(x3 <= 0.0 || x3 >= 0.33) {
-    nwarnings++;
-    if(nwarnings < nwarnlimit) 
-      G4cout << "G4MuMscModel::SampleScattering: ePre(MeV)= " << preKinEnergy/MeV
-	     << " ePost(MeV)= " << kinEnergy/MeV
-	     << " <x>= " << x1 << " sqrt(<x^2>)= " << sqrt(x2)
-	     << " x3= " << x3
-	     << G4endl;
-    x3 = std::min(1.0/y1,0.16666);
-  }
-  G4double x4 = 0.25*(3.0*x3 + sqrt(x3*(x3 + 8.0)))/(1.0 - x3);
-  */
-
-  G4double x  = G4UniformRand();
-  G4double z;
-  
-  //if(x < y1) z = x1*pow(x/y1,x4);
-  //else       z = 1.0 - y1*pow((1.0 - x)/x1,x4);
- 
-  z = -x1*log(x);
-
-  G4double cost = 1.0 - 2.0*z;
-  if(cost < -1.0) cost = -1.0;
-  else if(cost > 1.0) cost = 1.0;
-  G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
-
-  G4double phi  = twopi*G4UniformRand();
-
-  G4double dirx = sint*cos(phi);
-  G4double diry = sint*sin(phi);
-
-  //  G4cout << "G4MuMscModel::SampleSecondaries: tstep(mm)= " << truestep/mm
-  //	 << " lambdaeff= " << lambdaeff
-  //	 << " rms= " << rms << G4endl;
-
-  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
-  G4ThreeVector newDirection(dirx,diry,cost);
-  newDirection.rotateUz(oldDirection);
-  fParticleChange->ProposeMomentumDirection(newDirection);
-
-  if (latDisplasment && safety > tlimitminfix) {
-    G4double rms= sqrt(2.0*x1);
-    G4double rx = zPathLength*(0.5*dirx + invsqrt12*G4RandGauss::shoot(0.0,rms));
-    G4double ry = zPathLength*(0.5*diry + invsqrt12*G4RandGauss::shoot(0.0,rms));
-    G4double r  = sqrt(rx*rx + ry*ry);
-/*
-    G4cout << "G4MuMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
-	   << " sinTheta= " << sth << " r(mm)= " << r
-           << " trueStep(mm)= " << truestep 
-           << " geomStep(mm)= " << zPathLength
-           << G4endl;
-*/
-    
-    G4ThreeVector latDirection(rx,ry,0.0);
-    latDirection.rotateUz(oldDirection);
-
-    G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
-    G4double fac = 1.;
-    if(r >  safety) {
-      //  ******* so safety is computed at boundary too ************
-      G4double newsafety = safetyHelper->ComputeSafety(Position);
-      if(r > newsafety)
-	fac = newsafety/r ;
-    }  
-
-    if(fac > 0.) {
-      // compute new endpoint of the Step
-      G4ThreeVector newPosition = Position+fac*r*latDirection;
-
-      // definitely not on boundary
-      if(1. == fac) {
-	safetyHelper->ReLocateWithinVolume(newPosition);
-	    
-      } else {
-	// check safety after displacement
-	G4double postsafety = safetyHelper->ComputeSafety(newPosition);
-
-	// displacement to boundary
-	if(postsafety <= 0.0) {
-	  safetyHelper->Locate(newPosition, newDirection);
-
-	  // not on the boundary
-	} else { 
-	      safetyHelper->ReLocateWithinVolume(newPosition);
-	}
-      }
-      fParticleChange->ProposePosition(newPosition);
-    } 
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-void G4MuMscModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
-                                     const G4MaterialCutsCouple*,
-				     const G4DynamicParticle*,
-				     G4double,
-				     G4double)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
