Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4UrbanMscModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4UrbanMscModel.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4UrbanMscModel.cc (Version 11.2.1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // -------------------------------------------     26 // -------------------------------------------------------------------
 27 //                                                 27 //   
 28 // GEANT4 Class file                               28 // GEANT4 Class file
 29 //                                                 29 //    
 30 //                                                 30 //
 31 // File name:   G4UrbanMscModel                    31 // File name:   G4UrbanMscModel
 32 //                                                 32 //
 33 // Author:      Laszlo Urban                       33 // Author:      Laszlo Urban
 34 //                                                 34 //
 35 // Creation date: 19.02.2013                       35 // Creation date: 19.02.2013
 36 //                                                 36 //
 37 // Created from G4UrbanMscModel96                  37 // Created from G4UrbanMscModel96
 38 //                                                 38 //
 39 // New parametrization for theta0                  39 // New parametrization for theta0
 40 // Correction for very small step length           40 // Correction for very small step length
 41 //                                                 41 //
 42 // Class Description:                              42 // Class Description:
 43 //                                                 43 //
 44 // Implementation of the model of multiple sca     44 // Implementation of the model of multiple scattering based on
 45 // H.W.Lewis Phys Rev 78 (1950) 526 and others     45 // H.W.Lewis Phys Rev 78 (1950) 526 and others
 46                                                    46 
 47 // -------------------------------------------     47 // -------------------------------------------------------------------
 48 // In its present form the model can be  used      48 // In its present form the model can be  used for simulation 
 49 //   of the e-/e+ multiple scattering              49 //   of the e-/e+ multiple scattering
 50                                                    50 
 51                                                    51 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53 //....oooOO0OOooo........oooOO0OOooo........oo     53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 54                                                    54 
 55 #include "G4UrbanMscModel.hh"                      55 #include "G4UrbanMscModel.hh"
 56 #include "G4PhysicalConstants.hh"                  56 #include "G4PhysicalConstants.hh"
 57 #include "G4SystemOfUnits.hh"                      57 #include "G4SystemOfUnits.hh"
 58 #include "Randomize.hh"                            58 #include "Randomize.hh"
 59 #include "G4Positron.hh"                           59 #include "G4Positron.hh"
 60 #include "G4EmParameters.hh"                       60 #include "G4EmParameters.hh"
 61 #include "G4ParticleChangeForMSC.hh"               61 #include "G4ParticleChangeForMSC.hh"
 62 #include "G4ProductionCutsTable.hh"                62 #include "G4ProductionCutsTable.hh"
 63                                                    63 
 64 #include "G4Poisson.hh"                            64 #include "G4Poisson.hh"
 65 #include "G4Pow.hh"                                65 #include "G4Pow.hh"
 66 #include "G4Log.hh"                                66 #include "G4Log.hh"
 67 #include "G4Exp.hh"                                67 #include "G4Exp.hh"
 68 #include "G4AutoLock.hh"                           68 #include "G4AutoLock.hh"
 69                                                    69 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 71                                                    71 
 72 std::vector<G4UrbanMscModel::mscData*> G4Urban     72 std::vector<G4UrbanMscModel::mscData*> G4UrbanMscModel::msc;
 73                                                    73 
 74 namespace                                          74 namespace
 75 {                                                  75 {
 76   G4Mutex theUrbanMutex = G4MUTEX_INITIALIZER;     76   G4Mutex theUrbanMutex = G4MUTEX_INITIALIZER;
 77 }                                                  77 }
 78                                                    78 
 79 //....oooOO0OOooo........oooOO0OOooo........oo     79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 80                                                    80 
 81 G4UrbanMscModel::G4UrbanMscModel(const G4Strin     81 G4UrbanMscModel::G4UrbanMscModel(const G4String& nam)
 82   : G4VMscModel(nam)                               82   : G4VMscModel(nam)
 83 {                                                  83 {
 84   masslimite    = 0.6*CLHEP::MeV;                  84   masslimite    = 0.6*CLHEP::MeV;
 85   fr            = 0.02;                            85   fr            = 0.02;
 86   taubig        = 8.0;                             86   taubig        = 8.0;
 87   tausmall      = 1.e-16;                          87   tausmall      = 1.e-16;
 88   taulim        = 1.e-6;                           88   taulim        = 1.e-6;
 89   currentTau    = taulim;                          89   currentTau    = taulim;
 90   tlimitminfix  = 0.01*CLHEP::nm;                  90   tlimitminfix  = 0.01*CLHEP::nm;             
 91   tlimitminfix2 = 1.*CLHEP::nm;                    91   tlimitminfix2 = 1.*CLHEP::nm;
 92   stepmin       = tlimitminfix;                    92   stepmin       = tlimitminfix;
 93   smallstep     = 1.e10;                           93   smallstep     = 1.e10;
 94   currentRange  = 0.;                              94   currentRange  = 0.;
 95   rangeinit     = 0.;                              95   rangeinit     = 0.;
 96   tlimit        = 1.e10*CLHEP::mm;                 96   tlimit        = 1.e10*CLHEP::mm;
 97   tlimitmin     = 10.*tlimitminfix;                97   tlimitmin     = 10.*tlimitminfix;            
 98   tgeom         = 1.e50*CLHEP::mm;                 98   tgeom         = 1.e50*CLHEP::mm;
 99   geombig       = tgeom;                           99   geombig       = tgeom;
100   geommin       = 1.e-3*CLHEP::mm;                100   geommin       = 1.e-3*CLHEP::mm;
101   geomlimit     = geombig;                        101   geomlimit     = geombig;
102   presafety     = 0.;                             102   presafety     = 0.;
103                                                   103 
104   positron      = G4Positron::Positron();         104   positron      = G4Positron::Positron();
105   rndmEngineMod = G4Random::getTheEngine();       105   rndmEngineMod = G4Random::getTheEngine();
106                                                   106 
107   drr = 0.35;                                     107   drr = 0.35;
108   finalr = 10.*CLHEP::um;                         108   finalr = 10.*CLHEP::um;
109                                                   109 
110   tlow = 5.*CLHEP::keV;                           110   tlow = 5.*CLHEP::keV;
111   invmev = 1.0/CLHEP::MeV;                        111   invmev = 1.0/CLHEP::MeV;
112                                                   112 
113   skindepth = skin*stepmin;                       113   skindepth = skin*stepmin;
114                                                   114 
115   mass = CLHEP::proton_mass_c2;                   115   mass = CLHEP::proton_mass_c2;
116   charge = chargeSquare = 1.0;                    116   charge = chargeSquare = 1.0;
117   currentKinEnergy = currentRadLength = lambda    117   currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength 
118     = zPathLength = par1 = par2 = par3 = rndma    118     = zPathLength = par1 = par2 = par3 = rndmarray[0] = rndmarray[1] = 0;
119   currentLogKinEnergy = LOG_EKIN_MIN;             119   currentLogKinEnergy = LOG_EKIN_MIN;
120 }                                                 120 }
121                                                   121 
122 //....oooOO0OOooo........oooOO0OOooo........oo    122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123                                                   123 
124 G4UrbanMscModel::~G4UrbanMscModel()               124 G4UrbanMscModel::~G4UrbanMscModel()
125 {                                                 125 {
126   if(isFirstInstance) {                           126   if(isFirstInstance) {
127     for(auto const & ptr : msc) { delete ptr;  << 127     for(auto & ptr : msc) { delete ptr; }
128     msc.clear();                                  128     msc.clear();
129   }                                               129   } 
130 }                                                 130 }
131                                                   131 
132 //....oooOO0OOooo........oooOO0OOooo........oo    132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133                                                   133 
134 void G4UrbanMscModel::Initialise(const G4Parti    134 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p,
135                                  const G4DataV    135                                  const G4DataVector&)
136 {                                                 136 {
137   // set values of some data members              137   // set values of some data members
138   SetParticle(p);                                 138   SetParticle(p);
139   fParticleChange = GetParticleChangeForMSC(p)    139   fParticleChange = GetParticleChangeForMSC(p);
140   InitialiseParameters(p);                        140   InitialiseParameters(p);
141                                                   141 
142   latDisplasmentbackup = latDisplasment;          142   latDisplasmentbackup = latDisplasment;
143                                                   143 
144   // if model is locked parameters should be d    144   // if model is locked parameters should be defined via Set methods
145   if(!IsLocked()) {                               145   if(!IsLocked()) {
146     dispAlg96 = G4EmParameters::Instance()->La    146     dispAlg96 = G4EmParameters::Instance()->LateralDisplacementAlg96();
147     fPosiCorrection = G4EmParameters::Instance    147     fPosiCorrection = G4EmParameters::Instance()->MscPositronCorrection();
148   }                                               148   }
149                                                   149 
150   // initialise cache only once                   150   // initialise cache only once
151   if(0 == msc.size()) {                           151   if(0 == msc.size()) {
152     G4AutoLock l(&theUrbanMutex);                 152     G4AutoLock l(&theUrbanMutex);
153     if(0 == msc.size()) {                         153     if(0 == msc.size()) {
154       isFirstInstance = true;                     154       isFirstInstance = true;
155       msc.resize(1, nullptr);                     155       msc.resize(1, nullptr);
156     }                                             156     }
157     l.unlock();                                   157     l.unlock();
158   }                                               158   }
159   // initialise cache for each new run            159   // initialise cache for each new run
160   if(isFirstInstance) { InitialiseModelCache()    160   if(isFirstInstance) { InitialiseModelCache(); }
161                                                   161 
162   /*                                              162   /*
163   G4cout << "### G4UrbanMscModel::Initialise d    163   G4cout << "### G4UrbanMscModel::Initialise done for " 
164    << p->GetParticleName() << " type= " << ste    164    << p->GetParticleName() << " type= " << steppingAlgorithm << G4endl;
165   G4cout << "    RangeFact= " << facrange << "    165   G4cout << "    RangeFact= " << facrange << " GeomFact= " << facgeom
166    << " SafetyFact= " << facsafety << " Lambda    166    << " SafetyFact= " << facsafety << " LambdaLim= " << lambdalimit
167    << G4endl;                                     167    << G4endl;
168   */                                              168   */
169 }                                                 169 }
170                                                   170 
171 //....oooOO0OOooo........oooOO0OOooo........oo    171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172                                                   172 
173 G4double G4UrbanMscModel::ComputeCrossSectionP    173 G4double G4UrbanMscModel::ComputeCrossSectionPerAtom( 
174                              const G4ParticleD    174                              const G4ParticleDefinition* part,
175                                    G4double ki    175                                    G4double kinEnergy,
176                                    G4double at    176                                    G4double atomicNumber,G4double,
177                                    G4double, G    177                                    G4double, G4double)
178 {                                                 178 {
179   static const G4double epsmin = 1.e-4 , epsma    179   static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
180                                                   180 
181   static const G4double Zdat[15] = { 4.,  6.,     181   static const G4double Zdat[15] = { 4.,  6., 13., 20., 26., 29., 32., 38.,47.,
182                                      50., 56.,    182                                      50., 56., 64., 74., 79., 82. };
183                                                   183 
184   // corr. factors for e-/e+ lambda for T <= T    184   // corr. factors for e-/e+ lambda for T <= Tlim
185   static const G4double celectron[15][22] =       185   static const G4double celectron[15][22] =
186           {{1.125,1.072,1.051,1.047,1.047,1.05    186           {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
187             1.054,1.057,1.062,1.069,1.075,1.09    187             1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
188             1.112,1.108,1.100,1.093,1.089,1.08    188             1.112,1.108,1.100,1.093,1.089,1.087            },
189            {1.408,1.246,1.143,1.096,1.077,1.05    189            {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
190             1.052,1.053,1.058,1.065,1.072,1.08    190             1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
191             1.109,1.105,1.097,1.090,1.086,1.08    191             1.109,1.105,1.097,1.090,1.086,1.082            },
192            {2.833,2.268,1.861,1.612,1.486,1.30    192            {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
193             1.136,1.114,1.106,1.106,1.109,1.11    193             1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
194             1.131,1.124,1.113,1.104,1.099,1.09    194             1.131,1.124,1.113,1.104,1.099,1.098            },
195            {3.879,3.016,2.380,2.007,1.818,1.53    195            {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
196             1.190,1.133,1.107,1.099,1.098,1.10    196             1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
197             1.112,1.105,1.096,1.089,1.085,1.09    197             1.112,1.105,1.096,1.089,1.085,1.098            },
198            {6.937,4.330,2.886,2.256,1.987,1.62    198            {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
199             1.203,1.122,1.080,1.065,1.061,1.06    199             1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
200             1.073,1.070,1.064,1.059,1.056,1.05    200             1.073,1.070,1.064,1.059,1.056,1.056            },
201            {9.616,5.708,3.424,2.551,2.204,1.76    201            {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
202             1.256,1.155,1.099,1.077,1.070,1.06    202             1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
203             1.074,1.070,1.063,1.059,1.056,1.05    203             1.074,1.070,1.063,1.059,1.056,1.052            },
204            {11.72,6.364,3.811,2.806,2.401,1.88    204            {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
205             1.300,1.180,1.112,1.082,1.073,1.06    205             1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
206             1.068,1.064,1.059,1.054,1.051,1.05    206             1.068,1.064,1.059,1.054,1.051,1.050            },
207            {18.08,8.601,4.569,3.183,2.662,2.02    207            {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
208             1.339,1.195,1.108,1.068,1.053,1.04    208             1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
209             1.039,1.037,1.034,1.031,1.030,1.03    209             1.039,1.037,1.034,1.031,1.030,1.036            },
210            {18.22,10.48,5.333,3.713,3.115,2.36    210            {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
211             1.498,1.301,1.171,1.105,1.077,1.04    211             1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
212             1.031,1.028,1.024,1.022,1.021,1.02    212             1.031,1.028,1.024,1.022,1.021,1.024            },
213            {14.14,10.65,5.710,3.929,3.266,2.45    213            {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
214             1.528,1.319,1.178,1.106,1.075,1.04    214             1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
215             1.020,1.017,1.015,1.013,1.013,1.02    215             1.020,1.017,1.015,1.013,1.013,1.020            },
216            {14.11,11.73,6.312,4.240,3.478,2.56    216            {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
217             1.569,1.342,1.186,1.102,1.065,1.02    217             1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
218             0.995,0.993,0.993,0.993,0.993,1.01    218             0.995,0.993,0.993,0.993,0.993,1.011            },
219            {22.76,20.01,8.835,5.287,4.144,2.90    219            {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
220             1.677,1.410,1.224,1.121,1.073,1.01    220             1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
221             0.974,0.972,0.973,0.974,0.975,0.98    221             0.974,0.972,0.973,0.974,0.975,0.987            },
222            {50.77,40.85,14.13,7.184,5.284,3.43    222            {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
223             1.837,1.512,1.283,1.153,1.091,1.01    223             1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
224             0.950,0.947,0.949,0.952,0.954,0.96    224             0.950,0.947,0.949,0.952,0.954,0.963            },
225            {65.87,59.06,15.87,7.570,5.567,3.65    225            {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
226             1.939,1.579,1.325,1.178,1.108,1.01    226             1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
227             0.941,0.938,0.940,0.944,0.946,0.95    227             0.941,0.938,0.940,0.944,0.946,0.954            },
228            {55.60,47.34,15.92,7.810,5.755,3.76    228            {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
229             1.985,1.609,1.343,1.188,1.113,1.01    229             1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
230             0.933,0.930,0.933,0.936,0.939,0.94    230             0.933,0.930,0.933,0.936,0.939,0.949            }};
231                                                   231             
232   static const G4double cpositron[15][22] = {     232   static const G4double cpositron[15][22] = {
233            {2.589,2.044,1.658,1.446,1.347,1.21    233            {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
234             1.097,1.083,1.080,1.086,1.092,1.10    234             1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
235             1.131,1.126,1.117,1.108,1.103,1.10    235             1.131,1.126,1.117,1.108,1.103,1.100            },
236            {3.904,2.794,2.079,1.710,1.543,1.32    236            {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
237             1.122,1.096,1.089,1.092,1.098,1.11    237             1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
238             1.138,1.132,1.122,1.113,1.108,1.10    238             1.138,1.132,1.122,1.113,1.108,1.102            },
239            {7.970,6.080,4.442,3.398,2.872,2.12    239            {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
240             1.357,1.246,1.194,1.179,1.178,1.18    240             1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
241             1.203,1.190,1.173,1.159,1.151,1.14    241             1.203,1.190,1.173,1.159,1.151,1.145            },
242            {9.714,7.607,5.747,4.493,3.815,2.77    242            {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
243             1.553,1.353,1.253,1.219,1.211,1.21    243             1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
244             1.225,1.210,1.191,1.175,1.166,1.17    244             1.225,1.210,1.191,1.175,1.166,1.174            },
245            {17.97,12.95,8.628,6.065,4.849,3.22    245            {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
246             1.624,1.382,1.259,1.214,1.202,1.20    246             1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
247             1.217,1.203,1.184,1.169,1.160,1.15    247             1.217,1.203,1.184,1.169,1.160,1.151            },
248            {24.83,17.06,10.84,7.355,5.767,3.70    248            {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
249             1.759,1.465,1.311,1.252,1.234,1.22    249             1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
250             1.237,1.222,1.201,1.184,1.174,1.15    250             1.237,1.222,1.201,1.184,1.174,1.159            },
251            {23.26,17.15,11.52,8.049,6.375,4.11    251            {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
252             1.880,1.535,1.353,1.281,1.258,1.24    252             1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
253             1.252,1.234,1.212,1.194,1.183,1.17    253             1.252,1.234,1.212,1.194,1.183,1.170            },
254            {22.33,18.01,12.86,9.212,7.336,4.70    254            {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
255             2.015,1.602,1.385,1.297,1.268,1.25    255             2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
256             1.254,1.237,1.214,1.195,1.185,1.17    256             1.254,1.237,1.214,1.195,1.185,1.179            },
257            {33.91,24.13,15.71,10.80,8.507,5.46    257            {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
258             2.407,1.873,1.564,1.425,1.374,1.33    258             2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
259             1.312,1.288,1.258,1.235,1.221,1.20    259             1.312,1.288,1.258,1.235,1.221,1.205            },
260            {32.14,24.11,16.30,11.40,9.015,5.78    260            {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
261             2.490,1.925,1.596,1.447,1.391,1.34    261             2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
262             1.320,1.294,1.264,1.240,1.226,1.21    262             1.320,1.294,1.264,1.240,1.226,1.214            },
263            {29.51,24.07,17.19,12.28,9.766,6.23    263            {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
264             2.602,1.995,1.641,1.477,1.414,1.35    264             2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
265             1.328,1.302,1.270,1.245,1.231,1.23    265             1.328,1.302,1.270,1.245,1.231,1.233            },
266            {38.19,30.85,21.76,15.35,12.07,7.52    266            {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
267             2.926,2.188,1.763,1.563,1.484,1.40    267             2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
268             1.361,1.330,1.294,1.267,1.251,1.23    268             1.361,1.330,1.294,1.267,1.251,1.239            },
269            {49.71,39.80,27.96,19.63,15.36,9.40    269            {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
270             3.417,2.478,1.944,1.692,1.589,1.48    270             3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
271             1.409,1.372,1.330,1.298,1.280,1.25    271             1.409,1.372,1.330,1.298,1.280,1.258            },
272            {59.25,45.08,30.36,20.83,16.15,9.83    272            {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
273             3.641,2.648,2.064,1.779,1.661,1.53    273             3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
274             1.442,1.400,1.354,1.319,1.299,1.27    274             1.442,1.400,1.354,1.319,1.299,1.272            },
275            {56.38,44.29,30.50,21.18,16.51,10.1    275            {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
276             3.752,2.724,2.116,1.817,1.692,1.55    276             3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
277             1.456,1.412,1.364,1.328,1.307,1.28    277             1.456,1.412,1.364,1.328,1.307,1.282            }};
278                                                   278 
279   //data/corrections for T > Tlim                 279   //data/corrections for T > Tlim  
280                                                   280                                              
281   static const G4double hecorr[15] = {            281   static const G4double hecorr[15] = {
282     120.70, 117.50, 105.00, 92.92, 79.23,  74.    282     120.70, 117.50, 105.00, 92.92, 79.23,  74.510,  68.29,
283     57.39,  41.97,  36.14, 24.53, 10.21,  -7.8    283     57.39,  41.97,  36.14, 24.53, 10.21,  -7.855, -16.84,
284     -22.30};                                      284     -22.30};
285                                                   285 
286   G4double sigma;                                 286   G4double sigma;
287   SetParticle(part);                              287   SetParticle(part);
288                                                   288 
289   G4double Z23 = G4Pow::GetInstance()->Z23(G4l    289   G4double Z23 = G4Pow::GetInstance()->Z23(G4lrint(atomicNumber));
290                                                   290 
291   // correction if particle .ne. e-/e+            291   // correction if particle .ne. e-/e+
292   // compute equivalent kinetic energy            292   // compute equivalent kinetic energy
293   // lambda depends on p*beta ....                293   // lambda depends on p*beta ....
294                                                   294 
295   G4double eKineticEnergy = kinEnergy;            295   G4double eKineticEnergy = kinEnergy;
296                                                   296 
297   if(mass > CLHEP::electron_mass_c2)              297   if(mass > CLHEP::electron_mass_c2)
298   {                                               298   {
299      G4double TAU = kinEnergy/mass ;              299      G4double TAU = kinEnergy/mass ;
300      G4double c = mass*TAU*(TAU+2.)/(CLHEP::el    300      G4double c = mass*TAU*(TAU+2.)/(CLHEP::electron_mass_c2*(TAU+1.)) ;
301      G4double w = c-2.;                           301      G4double w = c-2.;
302      G4double tau = 0.5*(w+std::sqrt(w*w+4.*c)    302      G4double tau = 0.5*(w+std::sqrt(w*w+4.*c)) ;
303      eKineticEnergy = CLHEP::electron_mass_c2*    303      eKineticEnergy = CLHEP::electron_mass_c2*tau ;
304   }                                               304   }
305                                                   305 
306   G4double eTotalEnergy = eKineticEnergy + CLH    306   G4double eTotalEnergy = eKineticEnergy + CLHEP::electron_mass_c2 ;
307   G4double beta2 = eKineticEnergy*(eTotalEnerg    307   G4double beta2 = eKineticEnergy*(eTotalEnergy+CLHEP::electron_mass_c2)
308                                  /(eTotalEnerg    308                                  /(eTotalEnergy*eTotalEnergy);
309   G4double bg2   = eKineticEnergy*(eTotalEnerg    309   G4double bg2   = eKineticEnergy*(eTotalEnergy+CLHEP::electron_mass_c2)
310                                  /(CLHEP::elec    310                                  /(CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
311                                                   311 
312   static const G4double epsfactor = 2.*CLHEP::    312   static const G4double epsfactor = 2.*CLHEP::electron_mass_c2*
313     CLHEP::electron_mass_c2*CLHEP::Bohr_radius    313     CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
314     /(CLHEP::hbarc*CLHEP::hbarc);                 314     /(CLHEP::hbarc*CLHEP::hbarc);
315   G4double eps = epsfactor*bg2/Z23;               315   G4double eps = epsfactor*bg2/Z23;
316                                                   316 
317   if     (eps<epsmin)  sigma = 2.*eps*eps;        317   if     (eps<epsmin)  sigma = 2.*eps*eps;
318   else if(eps<epsmax)  sigma = G4Log(1.+2.*eps    318   else if(eps<epsmax)  sigma = G4Log(1.+2.*eps)-2.*eps/(1.+2.*eps);
319   else                 sigma = G4Log(2.*eps)-1    319   else                 sigma = G4Log(2.*eps)-1.+1./eps;
320                                                   320 
321   sigma *= chargeSquare*atomicNumber*atomicNum    321   sigma *= chargeSquare*atomicNumber*atomicNumber/(beta2*bg2);
322                                                   322 
323   // interpolate in AtomicNumber and beta2        323   // interpolate in AtomicNumber and beta2 
324   G4double c1,c2,cc1;                             324   G4double c1,c2,cc1;
325                                                   325 
326   // get bin number in Z                          326   // get bin number in Z
327   G4int iZ = 14;                                  327   G4int iZ = 14;
328   // Loop checking, 03-Aug-2015, Vladimir Ivan    328   // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
329   while ((iZ>=0)&&(Zdat[iZ]>=atomicNumber)) {     329   while ((iZ>=0)&&(Zdat[iZ]>=atomicNumber)) { --iZ; }
330                                                   330 
331   iZ = std::min(std::max(iZ, 0), 13);             331   iZ = std::min(std::max(iZ, 0), 13);
332                                                   332 
333   G4double ZZ1 = Zdat[iZ];                        333   G4double ZZ1 = Zdat[iZ];
334   G4double ZZ2 = Zdat[iZ+1];                      334   G4double ZZ2 = Zdat[iZ+1];
335   G4double ratZ = (atomicNumber-ZZ1)*(atomicNu    335   G4double ratZ = (atomicNumber-ZZ1)*(atomicNumber+ZZ1)/
336                   ((ZZ2-ZZ1)*(ZZ2+ZZ1));          336                   ((ZZ2-ZZ1)*(ZZ2+ZZ1));
337                                                   337 
338   static const G4double Tlim = 10.*CLHEP::MeV;    338   static const G4double Tlim = 10.*CLHEP::MeV;
339   static const G4double sigmafactor =             339   static const G4double sigmafactor =
340     CLHEP::twopi*CLHEP::classic_electr_radius*    340     CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
341   static const G4double beta2lim = Tlim*(Tlim+    341   static const G4double beta2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
342     ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHE    342     ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHEP::electron_mass_c2));
343   static const G4double bg2lim   = Tlim*(Tlim+    343   static const G4double bg2lim   = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
344     (CLHEP::electron_mass_c2*CLHEP::electron_m    344     (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
345                                                   345 
346   static const G4double sig0[15] = {              346   static const G4double sig0[15] = {
347     0.2672*CLHEP::barn,  0.5922*CLHEP::barn,      347     0.2672*CLHEP::barn,  0.5922*CLHEP::barn,  2.653*CLHEP::barn, 6.235*CLHEP::barn,
348     11.69*CLHEP::barn  , 13.24*CLHEP::barn  ,     348     11.69*CLHEP::barn  , 13.24*CLHEP::barn  , 16.12*CLHEP::barn, 23.00*CLHEP::barn,
349     35.13*CLHEP::barn  , 39.95*CLHEP::barn  ,     349     35.13*CLHEP::barn  , 39.95*CLHEP::barn  , 50.85*CLHEP::barn, 67.19*CLHEP::barn,
350     91.15*CLHEP::barn  , 104.4*CLHEP::barn  ,     350     91.15*CLHEP::barn  , 104.4*CLHEP::barn  , 113.1*CLHEP::barn};
351                                                   351 
352   static const G4double Tdat[22] = {              352   static const G4double Tdat[22] = { 
353     100*CLHEP::eV,  200*CLHEP::eV,  400*CLHEP:    353     100*CLHEP::eV,  200*CLHEP::eV,  400*CLHEP::eV,  700*CLHEP::eV,
354     1*CLHEP::keV,   2*CLHEP::keV,   4*CLHEP::k    354     1*CLHEP::keV,   2*CLHEP::keV,   4*CLHEP::keV,   7*CLHEP::keV,
355     10*CLHEP::keV,  20*CLHEP::keV,  40*CLHEP::    355     10*CLHEP::keV,  20*CLHEP::keV,  40*CLHEP::keV,  70*CLHEP::keV,
356     100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP:    356     100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP::keV, 700*CLHEP::keV,
357     1*CLHEP::MeV,   2*CLHEP::MeV,   4*CLHEP::M    357     1*CLHEP::MeV,   2*CLHEP::MeV,   4*CLHEP::MeV,   7*CLHEP::MeV,
358     10*CLHEP::MeV,  20*CLHEP::MeV};               358     10*CLHEP::MeV,  20*CLHEP::MeV};
359                                                   359 
360   if(eKineticEnergy <= Tlim)                      360   if(eKineticEnergy <= Tlim) 
361   {                                               361   {
362     // get bin number in T (beta2)                362     // get bin number in T (beta2)
363     G4int iT = 21;                                363     G4int iT = 21;
364     // Loop checking, 03-Aug-2015, Vladimir Iv    364     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
365     while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)    365     while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
366                                                   366 
367     iT = std::min(std::max(iT, 0), 20);           367     iT = std::min(std::max(iT, 0), 20);
368                                                   368 
369     //  calculate betasquare values               369     //  calculate betasquare values
370     G4double T = Tdat[iT];                        370     G4double T = Tdat[iT];
371     G4double E = T + CLHEP::electron_mass_c2;     371     G4double E = T + CLHEP::electron_mass_c2;
372     G4double b2small = T*(E+CLHEP::electron_ma    372     G4double b2small = T*(E+CLHEP::electron_mass_c2)/(E*E);
373                                                   373 
374     T = Tdat[iT+1];                               374     T = Tdat[iT+1];
375     E = T + CLHEP::electron_mass_c2;              375     E = T + CLHEP::electron_mass_c2;
376     G4double b2big = T*(E+CLHEP::electron_mass    376     G4double b2big = T*(E+CLHEP::electron_mass_c2)/(E*E);
377     G4double ratb2 = (beta2-b2small)/(b2big-b2    377     G4double ratb2 = (beta2-b2small)/(b2big-b2small);
378                                                   378 
379     if (charge < 0.)                              379     if (charge < 0.)
380     {                                             380     {
381        c1 = celectron[iZ][iT];                    381        c1 = celectron[iZ][iT];
382        c2 = celectron[iZ+1][iT];                  382        c2 = celectron[iZ+1][iT];
383        cc1 = c1+ratZ*(c2-c1);                     383        cc1 = c1+ratZ*(c2-c1);
384                                                   384 
385        c1 = celectron[iZ][iT+1];                  385        c1 = celectron[iZ][iT+1];
386        c2 = celectron[iZ+1][iT+1];                386        c2 = celectron[iZ+1][iT+1];
387     }                                             387     }
388     else                                          388     else              
389     {                                             389     {
390        c1 = cpositron[iZ][iT];                    390        c1 = cpositron[iZ][iT];
391        c2 = cpositron[iZ+1][iT];                  391        c2 = cpositron[iZ+1][iT];
392        cc1 = c1+ratZ*(c2-c1);                     392        cc1 = c1+ratZ*(c2-c1);
393                                                   393 
394        c1 = cpositron[iZ][iT+1];                  394        c1 = cpositron[iZ][iT+1];
395        c2 = cpositron[iZ+1][iT+1];                395        c2 = cpositron[iZ+1][iT+1];
396     }                                             396     }
397     G4double cc2 = c1+ratZ*(c2-c1);               397     G4double cc2 = c1+ratZ*(c2-c1);
398     sigma *= sigmafactor/(cc1+ratb2*(cc2-cc1))    398     sigma *= sigmafactor/(cc1+ratb2*(cc2-cc1));
399   }                                               399   }
400   else                                            400   else
401   {                                               401   {
402     c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2    402     c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
403     c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(b    403     c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
404     if((atomicNumber >= ZZ1) && (atomicNumber     404     if((atomicNumber >= ZZ1) && (atomicNumber <= ZZ2))
405       sigma = c1+ratZ*(c2-c1) ;                   405       sigma = c1+ratZ*(c2-c1) ;
406     else if(atomicNumber < ZZ1)                   406     else if(atomicNumber < ZZ1)
407       sigma = atomicNumber*atomicNumber*c1/(ZZ    407       sigma = atomicNumber*atomicNumber*c1/(ZZ1*ZZ1);
408     else if(atomicNumber > ZZ2)                   408     else if(atomicNumber > ZZ2)
409       sigma = atomicNumber*atomicNumber*c2/(ZZ    409       sigma = atomicNumber*atomicNumber*c2/(ZZ2*ZZ2);
410   }                                               410   }
411   // low energy correction based on theory        411   // low energy correction based on theory 
412   sigma *= (1.+0.30/(1.+std::sqrt(1000.*eKinet    412   sigma *= (1.+0.30/(1.+std::sqrt(1000.*eKineticEnergy)));  
413                                                   413 
414   return sigma;                                   414   return sigma;
415 }                                                 415 }
416                                                   416 
417 //....oooOO0OOooo........oooOO0OOooo........oo    417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418                                                   418 
419 void G4UrbanMscModel::StartTracking(G4Track* t    419 void G4UrbanMscModel::StartTracking(G4Track* track)
420 {                                                 420 {
421   SetParticle(track->GetDynamicParticle()->Get    421   SetParticle(track->GetDynamicParticle()->GetDefinition());
422   firstStep = true;                               422   firstStep = true; 
423   insideskin = false;                             423   insideskin = false;
424   fr = facrange;                                  424   fr = facrange;
425   tlimit = tgeom = rangeinit = geombig;           425   tlimit = tgeom = rangeinit = geombig;
426   smallstep     = 1.e10;                          426   smallstep     = 1.e10;
427   stepmin       = tlimitminfix;                   427   stepmin       = tlimitminfix;
428   tlimitmin     = 10.*tlimitminfix;               428   tlimitmin     = 10.*tlimitminfix;            
429   rndmEngineMod = G4Random::getTheEngine();       429   rndmEngineMod = G4Random::getTheEngine();
430 }                                                 430 }
431                                                   431 
432 //....oooOO0OOooo........oooOO0OOooo........oo    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433                                                   433 
434 G4double G4UrbanMscModel::ComputeTruePathLengt    434 G4double G4UrbanMscModel::ComputeTruePathLengthLimit(
435                              const G4Track& tr    435                              const G4Track& track,
436                              G4double& current    436                              G4double& currentMinimalStep)
437 {                                                 437 {
438   tPathLength = currentMinimalStep;               438   tPathLength = currentMinimalStep;
439   const G4DynamicParticle* dp = track.GetDynam    439   const G4DynamicParticle* dp = track.GetDynamicParticle();
440                                                   440   
441   G4StepPoint* sp = track.GetStep()->GetPreSte    441   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
442   G4StepStatus stepStatus = sp->GetStepStatus(    442   G4StepStatus stepStatus = sp->GetStepStatus();
443   couple = track.GetMaterialCutsCouple();         443   couple = track.GetMaterialCutsCouple();
444   SetCurrentCouple(couple);                       444   SetCurrentCouple(couple); 
445   idx = couple->GetIndex();                       445   idx = couple->GetIndex();
446   currentKinEnergy = dp->GetKineticEnergy();      446   currentKinEnergy = dp->GetKineticEnergy();
447   currentLogKinEnergy = dp->GetLogKineticEnerg    447   currentLogKinEnergy = dp->GetLogKineticEnergy();
448   currentRange = GetRange(particle,currentKinE    448   currentRange = GetRange(particle,currentKinEnergy,couple,currentLogKinEnergy);
449   lambda0 = GetTransportMeanFreePath(particle,    449   lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy,
450                                                   450                                               currentLogKinEnergy);
451   tPathLength = std::min(tPathLength,currentRa    451   tPathLength = std::min(tPathLength,currentRange);
452   /*                                              452   /*  
453   G4cout << "G4Urban::StepLimit tPathLength= "    453   G4cout << "G4Urban::StepLimit tPathLength= " << tPathLength 
454   << " range= " <<currentRange<< " lambda= "<<    454   << " range= " <<currentRange<< " lambda= "<<lambda0
455             <<G4endl;                             455             <<G4endl;
456   */                                              456   */
457   // extreme small step                           457   // extreme small step
458   if(tPathLength < tlimitminfix) {                458   if(tPathLength < tlimitminfix) {
459     latDisplasment = false;                       459     latDisplasment = false;
460     return ConvertTrueToGeom(tPathLength, curr    460     return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
461   }                                               461   }
462                                                   462 
463   presafety = (stepStatus == fGeomBoundary) ?     463   presafety = (stepStatus == fGeomBoundary) ? sp->GetSafety()
464               : ComputeSafety(sp->GetPosition(    464               : ComputeSafety(sp->GetPosition(), tPathLength);
465                                                   465   
466   // stop here if small step or range is less     466   // stop here if small step or range is less than safety
467   if((tPathLength == currentRange && tPathLeng    467   if((tPathLength == currentRange && tPathLength < presafety) ||
468      tPathLength < tlimitminfix) {                468      tPathLength < tlimitminfix) { 
469     latDisplasment = false;                       469     latDisplasment = false;   
470     return ConvertTrueToGeom(tPathLength, curr    470     return ConvertTrueToGeom(tPathLength, currentMinimalStep); 
471   }                                               471   }
472                                                   472 
473   // upper limit for the straight line distanc    473   // upper limit for the straight line distance the particle can travel
474   // for electrons and positrons                  474   // for electrons and positrons
475   G4double distance = (mass < masslimite)         475   G4double distance = (mass < masslimite) 
476     ? currentRange*msc[idx]->doverra              476     ? currentRange*msc[idx]->doverra
477     // for muons, hadrons                         477     // for muons, hadrons
478     : currentRange*msc[idx]->doverrb;             478     : currentRange*msc[idx]->doverrb;
479                                                   479 
480   /*                                              480   /*  
481   G4cout << "G4Urban::StepLimit tPathLength= "    481   G4cout << "G4Urban::StepLimit tPathLength= " 
482             <<tPathLength<<" safety= " << pres    482             <<tPathLength<<" safety= " << presafety
483           << " range= " <<currentRange<< " lam    483           << " range= " <<currentRange<< " lambda= "<<lambda0
484             << " Alg: " << steppingAlgorithm <    484             << " Alg: " << steppingAlgorithm <<G4endl;
485   */                                              485   */
486   // far from geometry boundary                   486   // far from geometry boundary
487   if(distance < presafety)                        487   if(distance < presafety)
488     {                                             488     {
489       latDisplasment = false;                     489       latDisplasment = false;   
490       return ConvertTrueToGeom(tPathLength, cu    490       return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
491     }                                             491     }
492                                                   492 
493   latDisplasment = latDisplasmentbackup;          493   latDisplasment = latDisplasmentbackup;
494   // -----------------------------------------    494   // ----------------------------------------------------------------
495   // distance to boundary                         495   // distance to boundary 
496   if (steppingAlgorithm == fUseDistanceToBound    496   if (steppingAlgorithm == fUseDistanceToBoundary)
497     {                                             497     {
498       //compute geomlimit and presafety           498       //compute geomlimit and presafety 
499       geomlimit = ComputeGeomLimit(track, pres    499       geomlimit = ComputeGeomLimit(track, presafety, currentRange);
500       /*                                          500       /*
501         G4cout << "G4Urban::Distance to bounda    501         G4cout << "G4Urban::Distance to boundary geomlimit= "
502             <<geomlimit<<" safety= " << presaf    502             <<geomlimit<<" safety= " << presafety<<G4endl;
503       */                                          503       */
504                                                   504 
505       smallstep += 1.;                            505       smallstep += 1.;
506       insideskin = false;                         506       insideskin = false;
507       tgeom = geombig;                         << 
508                                                   507 
509       // initialisation at first step and at t << 508       // initialisation at firs step and at the boundary
510       if(firstStep || (stepStatus == fGeomBoun    509       if(firstStep || (stepStatus == fGeomBoundary))
511         {                                         510         {
512           rangeinit = currentRange;               511           rangeinit = currentRange;
513           if(!firstStep) { smallstep = 1.; }      512           if(!firstStep) { smallstep = 1.; }
514                                                   513 
515           //stepmin ~ lambda_elastic              514           //stepmin ~ lambda_elastic
516           stepmin = ComputeStepmin();             515           stepmin = ComputeStepmin();
517           skindepth = skin*stepmin;               516           skindepth = skin*stepmin;
518           tlimitmin = ComputeTlimitmin();         517           tlimitmin = ComputeTlimitmin();
519         /*                                        518         /* 
520           G4cout << "rangeinit= " << rangeinit    519           G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
521                  << " tlimitmin= " << tlimitmi    520                  << " tlimitmin= " << tlimitmin << " geomlimit= " 
522                  << geomlimit <<G4endl;           521                  << geomlimit <<G4endl;
523         */                                        522         */
524         }                                      << 523           // constraint from the geometry
525       // constraint from the geometry          << 524 
526       if((geomlimit < geombig) && (geomlimit > << 525           if((geomlimit < geombig) && (geomlimit > geommin))
527         {                                      << 526             {
528           // geomlimit is a geometrical step l << 527               // geomlimit is a geometrical step length
529           // transform it to true path length  << 528               // transform it to true path length (estimation)
530           if(lambda0 > geomlimit) {            << 529               if(lambda0 > geomlimit) {
531             geomlimit = -lambda0*G4Log(1.-geom << 530                 geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin;
532           }                                    << 531         }
533           tgeom = (stepStatus == fGeomBoundary << 532               tgeom = (stepStatus == fGeomBoundary)
534       : facrange*rangeinit + stepmin;          << 533                 ? geomlimit/facgeom : 2.*geomlimit/facgeom;
                                                   >> 534             }
                                                   >> 535           else
                                                   >> 536       {
                                                   >> 537         tgeom = geombig;
                                                   >> 538       }
535         }                                         539         }
536                                                   540 
537       //step limit                                541       //step limit 
538       tlimit = (currentRange > presafety) ?       542       tlimit = (currentRange > presafety) ?
539         std::max(facrange*rangeinit, facsafety    543         std::max(facrange*rangeinit, facsafety*presafety) : currentRange;
540                                                   544 
541       //lower limit for tlimit                    545       //lower limit for tlimit
542       tlimit = std::min(std::max(tlimit,tlimit    546       tlimit = std::min(std::max(tlimit,tlimitmin), tgeom);
543       /*                                          547       /*
544       G4cout << "tgeom= " << tgeom << " geomli    548       G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
545             << " tlimit= " << tlimit << " pres    549             << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
546       */                                          550       */
547       // shortcut                                 551       // shortcut
548       if((tPathLength < tlimit) && (tPathLengt    552       if((tPathLength < tlimit) && (tPathLength < presafety) &&
549          (smallstep > skin) && (tPathLength <     553          (smallstep > skin) && (tPathLength < geomlimit-0.999*skindepth))
550       {                                           554       {
551         return ConvertTrueToGeom(tPathLength,     555         return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
552       }                                           556       }
553                                                   557 
554       // step reduction near to boundary          558       // step reduction near to boundary
555       if(smallstep <= skin)                       559       if(smallstep <= skin)
556         {                                         560         {
557           tlimit = stepmin;                       561           tlimit = stepmin;
558           insideskin = true;                      562           insideskin = true;
559         }                                         563         }
560       else if(geomlimit < geombig)                564       else if(geomlimit < geombig)
561         {                                         565         {
562           if(geomlimit > skindepth)               566           if(geomlimit > skindepth)
563             {                                     567             {
564               tlimit = std::min(tlimit, geomli    568               tlimit = std::min(tlimit, geomlimit-0.999*skindepth);
565             }                                     569             }
566           else                                    570           else
567             {                                     571             {
568               insideskin = true;                  572               insideskin = true;
569               tlimit = std::min(tlimit, stepmi    573               tlimit = std::min(tlimit, stepmin);
570             }                                     574             }
571         }                                         575         }
572                                                   576 
573       tlimit = std::max(tlimit, stepmin);         577       tlimit = std::max(tlimit, stepmin); 
574                                                   578 
575       // randomise if not 'small' step and ste    579       // randomise if not 'small' step and step determined by msc
576       tPathLength = (tlimit < tPathLength && s    580       tPathLength = (tlimit < tPathLength && smallstep > skin && !insideskin) 
577         ? std::min(tPathLength, Randomizetlimi    581         ? std::min(tPathLength, Randomizetlimit()) 
578   : std::min(tPathLength, tlimit);                582   : std::min(tPathLength, tlimit);
579     }                                             583     }
580   // -----------------------------------------    584   // ----------------------------------------------------------------
581   // for simulation with or without magnetic f    585   // for simulation with or without magnetic field 
582   // there no small step/single scattering at     586   // there no small step/single scattering at boundaries
583   else if(steppingAlgorithm == fUseSafety)        587   else if(steppingAlgorithm == fUseSafety)
584     {                                             588     {
585       if(firstStep || (stepStatus == fGeomBoun    589       if(firstStep || (stepStatus == fGeomBoundary)) {
586         rangeinit = currentRange;                 590         rangeinit = currentRange;
587         fr = facrange;                            591         fr = facrange;
588         // stepping for e+/e- only (not for mu    592         // stepping for e+/e- only (not for muons,hadrons)
589         if(mass < masslimite)                     593         if(mass < masslimite) 
590           {                                       594           {
591             rangeinit = std::max(rangeinit, la    595             rangeinit = std::max(rangeinit, lambda0);
592             if(lambda0 > lambdalimit) {           596             if(lambda0 > lambdalimit) {
593               fr *= (0.75+0.25*lambda0/lambdal    597               fr *= (0.75+0.25*lambda0/lambdalimit);
594             }                                     598             }
595           }                                       599           }
596         //lower limit for tlimit                  600         //lower limit for tlimit
597         stepmin = ComputeStepmin();               601         stepmin = ComputeStepmin();
598         tlimitmin = ComputeTlimitmin();           602         tlimitmin = ComputeTlimitmin();
599       }                                           603       }
600                                                   604 
601       //step limit                                605       //step limit
602       tlimit = (currentRange > presafety) ?       606       tlimit = (currentRange > presafety) ?
603         std::max(fr*rangeinit, facsafety*presa    607         std::max(fr*rangeinit, facsafety*presafety) : currentRange;
604                                                   608   
605       //lower limit for tlimit                    609       //lower limit for tlimit
606       tlimit = std::max(tlimit, tlimitmin);       610       tlimit = std::max(tlimit, tlimitmin); 
607                                                   611      
608       // randomise if step determined by msc      612       // randomise if step determined by msc
609       tPathLength = (tlimit < tPathLength) ?      613       tPathLength = (tlimit < tPathLength) ?
610         std::min(tPathLength, Randomizetlimit(    614         std::min(tPathLength, Randomizetlimit()) : tPathLength; 
611     }                                             615     }
612   // -----------------------------------------    616   // ----------------------------------------------------------------
613   // for simulation with or without magnetic f    617   // for simulation with or without magnetic field 
614   // there is small step/single scattering at     618   // there is small step/single scattering at boundaries
615   else if(steppingAlgorithm == fUseSafetyPlus)    619   else if(steppingAlgorithm == fUseSafetyPlus)
616     {                                             620     {
617       if(firstStep || (stepStatus == fGeomBoun    621       if(firstStep || (stepStatus == fGeomBoundary)) {
618         rangeinit = currentRange;                 622         rangeinit = currentRange;
619         fr = facrange;                            623         fr = facrange;
620         if(mass < masslimite)                     624         if(mass < masslimite)
621           {                                       625           {
622             if(lambda0 > lambdalimit) {           626             if(lambda0 > lambdalimit) {
623               fr *= (0.84+0.16*lambda0/lambdal    627               fr *= (0.84+0.16*lambda0/lambdalimit);
624             }                                     628             }
625           }                                       629           }
626         //lower limit for tlimit                  630         //lower limit for tlimit
627         stepmin = ComputeStepmin();               631         stepmin = ComputeStepmin();
628         tlimitmin = ComputeTlimitmin();           632         tlimitmin = ComputeTlimitmin();
629       }                                           633       }
630       //step limit                                634       //step limit
631       tlimit = (currentRange > presafety) ?       635       tlimit = (currentRange > presafety) ?
632   std::max(fr*rangeinit, facsafety*presafety)     636   std::max(fr*rangeinit, facsafety*presafety) : currentRange;
633                                                   637 
634       //lower limit for tlimit                    638       //lower limit for tlimit
635       tlimit = std::max(tlimit, tlimitmin);       639       tlimit = std::max(tlimit, tlimitmin);
636                                                   640 
637       // condition for tPathLength from drr an    641       // condition for tPathLength from drr and finalr
638       if(currentRange > finalr) {                 642       if(currentRange > finalr) {
639         G4double tmax = drr*currentRange+         643         G4double tmax = drr*currentRange+
640                         finalr*(1.-drr)*(2.-fi    644                         finalr*(1.-drr)*(2.-finalr/currentRange);
641         tPathLength = std::min(tPathLength,tma    645         tPathLength = std::min(tPathLength,tmax); 
642       }                                           646       }
643                                                   647 
644       // randomise if step determined by msc      648       // randomise if step determined by msc
645       tPathLength = (tlimit < tPathLength) ?      649       tPathLength = (tlimit < tPathLength) ?
646         std::min(tPathLength, Randomizetlimit(    650         std::min(tPathLength, Randomizetlimit()) : tPathLength; 
647     }                                             651     }
648                                                   652 
649   // -----------------------------------------    653   // ----------------------------------------------------------------
650   // simple step limitation                       654   // simple step limitation
651   else                                            655   else
652     {                                             656     {
653       if (stepStatus == fGeomBoundary)            657       if (stepStatus == fGeomBoundary)
654         {                                         658         {
655           tlimit = (currentRange > lambda0)       659           tlimit = (currentRange > lambda0) 
656       ? facrange*currentRange : facrange*lambd    660       ? facrange*currentRange : facrange*lambda0;
657           tlimit = std::max(tlimit, tlimitmin)    661           tlimit = std::max(tlimit, tlimitmin);
658         }                                         662         }
659       // randomise if step determined by msc      663       // randomise if step determined by msc
660       tPathLength = (tlimit < tPathLength) ?      664       tPathLength = (tlimit < tPathLength) ?
661         std::min(tPathLength, Randomizetlimit(    665         std::min(tPathLength, Randomizetlimit()) : tPathLength; 
662     }                                             666     }
663                                                   667 
664   // -----------------------------------------    668   // ----------------------------------------------------------------
665   firstStep = false;                              669   firstStep = false; 
666   return ConvertTrueToGeom(tPathLength, curren    670   return ConvertTrueToGeom(tPathLength, currentMinimalStep);
667 }                                                 671 }
668                                                   672 
669 //....oooOO0OOooo........oooOO0OOooo........oo    673 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
670                                                   674 
671 G4double G4UrbanMscModel::ComputeGeomPathLengt    675 G4double G4UrbanMscModel::ComputeGeomPathLength(G4double)
672 {                                                 676 {
673   lambdaeff = lambda0;                            677   lambdaeff = lambda0;
674   par1 = -1. ;                                    678   par1 = -1. ;  
675   par2 = par3 = 0. ;                              679   par2 = par3 = 0. ;  
676                                                   680 
677   // this correction needed to run MSC with eI    681   // this correction needed to run MSC with eIoni and eBrem inactivated
678   // and makes no harm for a normal run           682   // and makes no harm for a normal run
679   tPathLength = std::min(tPathLength,currentRa    683   tPathLength = std::min(tPathLength,currentRange); 
680                                                   684 
681   //  do the true -> geom transformation          685   //  do the true -> geom transformation
682   zPathLength = tPathLength;                      686   zPathLength = tPathLength;
683                                                   687 
684   // z = t for very small tPathLength             688   // z = t for very small tPathLength
685   if(tPathLength < tlimitminfix2) return zPath    689   if(tPathLength < tlimitminfix2) return zPathLength;
686                                                   690 
687   /*                                              691   /*
688   G4cout << "ComputeGeomPathLength: tpl= " <<     692   G4cout << "ComputeGeomPathLength: tpl= " <<  tPathLength
689          << " R= " << currentRange << " L0= "     693          << " R= " << currentRange << " L0= " << lambda0
690          << " E= " << currentKinEnergy << "  "    694          << " E= " << currentKinEnergy << "  " 
691          << particle->GetParticleName() << G4e    695          << particle->GetParticleName() << G4endl;
692   */                                              696   */
693   G4double tau = tPathLength/lambda0 ;            697   G4double tau = tPathLength/lambda0 ;
694                                                   698 
695   if ((tau <= tausmall) || insideskin) {          699   if ((tau <= tausmall) || insideskin) {
696     zPathLength = std::min(tPathLength, lambda    700     zPathLength = std::min(tPathLength, lambda0); 
697                                                   701 
698   } else if (tPathLength < currentRange*dtrl)     702   } else if (tPathLength < currentRange*dtrl) {
699     zPathLength = (tau < taulim) ? tPathLength    703     zPathLength = (tau < taulim) ? tPathLength*(1.-0.5*tau)
700       : lambda0*(1.-G4Exp(-tau));                 704       : lambda0*(1.-G4Exp(-tau));
701                                                   705 
702   } else if(currentKinEnergy < mass || tPathLe    706   } else if(currentKinEnergy < mass || tPathLength == currentRange)  {
703     par1 = 1./currentRange;                       707     par1 = 1./currentRange;
704     par2 = currentRange/lambda0;                  708     par2 = currentRange/lambda0;
705     par3 = 1.+par2;                               709     par3 = 1.+par2;
706     if(tPathLength < currentRange) {              710     if(tPathLength < currentRange) {
707       zPathLength =                               711       zPathLength = 
708         (1.-G4Exp(par3*G4Log(1.-tPathLength/cu    712         (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3);
709     } else {                                      713     } else {
710       zPathLength = 1./(par1*par3);               714       zPathLength = 1./(par1*par3);
711     }                                             715     }
712                                                   716 
713   } else {                                        717   } else {
714     G4double rfin = std::max(currentRange-tPat    718     G4double rfin = std::max(currentRange-tPathLength, 0.01*currentRange);
715     G4double T1 = GetEnergy(particle,rfin,coup    719     G4double T1 = GetEnergy(particle,rfin,couple);
716     G4double lambda1 = GetTransportMeanFreePat    720     G4double lambda1 = GetTransportMeanFreePath(particle,T1);
717                                                   721 
718     par1 = (lambda0-lambda1)/(lambda0*tPathLen    722     par1 = (lambda0-lambda1)/(lambda0*tPathLength);
719     //G4cout << "par1= " << par1 << " L1= " <<    723     //G4cout << "par1= " << par1 << " L1= " << lambda1 << G4endl;
720     par2 = 1./(par1*lambda0);                     724     par2 = 1./(par1*lambda0);
721     par3 = 1.+par2;                               725     par3 = 1.+par2;
722     zPathLength = (1.-G4Exp(par3*G4Log(lambda1    726     zPathLength = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3);
723   }                                               727   }
724                                                   728 
725   zPathLength = std::min(zPathLength, lambda0)    729   zPathLength = std::min(zPathLength, lambda0);
726   //G4cout<< "zPathLength= "<< zPathLength<< "    730   //G4cout<< "zPathLength= "<< zPathLength<< " L0= " << lambda0 << G4endl;
727   return zPathLength;                             731   return zPathLength;
728 }                                                 732 }
729                                                   733 
730 //....oooOO0OOooo........oooOO0OOooo........oo    734 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
731                                                   735 
732 G4double G4UrbanMscModel::ComputeTrueStepLengt    736 G4double G4UrbanMscModel::ComputeTrueStepLength(G4double geomStepLength)
733 {                                                 737 {
734   // step defined other than transportation       738   // step defined other than transportation 
735   if(geomStepLength == zPathLength) {             739   if(geomStepLength == zPathLength) { 
736     //G4cout << "Urban::ComputeTrueLength: tPa    740     //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength 
737     //           << " step= " << geomStepLengt    741     //           << " step= " << geomStepLength << " *** " << G4endl;
738     return tPathLength;                           742     return tPathLength; 
739   }                                               743   }
740                                                   744 
741   zPathLength = geomStepLength;                   745   zPathLength = geomStepLength;
742                                                   746 
743   // t = z for very small step                    747   // t = z for very small step
744   if(geomStepLength < tlimitminfix2) {            748   if(geomStepLength < tlimitminfix2) { 
745     tPathLength = geomStepLength;                 749     tPathLength = geomStepLength; 
746                                                   750   
747   // recalculation                                751   // recalculation
748   } else {                                        752   } else {
749                                                   753 
750     G4double tlength = geomStepLength;            754     G4double tlength = geomStepLength;
751     if((geomStepLength > lambda0*tausmall) &&     755     if((geomStepLength > lambda0*tausmall) && !insideskin) {
752                                                   756 
753       if(par1 <  0.) {                            757       if(par1 <  0.) {
754         tlength = -lambda0*G4Log(1.-geomStepLe    758         tlength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
755       } else {                                    759       } else {
756         const G4double par4 = par1*par3;          760         const G4double par4 = par1*par3;
757         if(par4*geomStepLength < 1.) {            761         if(par4*geomStepLength < 1.) {
758           tlength = (1.-G4Exp(G4Log(1.-par4*ge    762           tlength = (1.-G4Exp(G4Log(1.-par4*geomStepLength)/par3))/par1;
759         } else {                                  763         } else {
760           tlength = currentRange;                 764           tlength = currentRange;
761         }                                         765         }
762       }                                           766       }
763                                                   767 
764       if(tlength < geomStepLength)   { tlength    768       if(tlength < geomStepLength)   { tlength = geomStepLength; }
765       else if(tlength > tPathLength) { tlength    769       else if(tlength > tPathLength) { tlength = tPathLength; }
766     }                                             770     }  
767     tPathLength = tlength;                        771     tPathLength = tlength; 
768   }                                               772   }
769   //G4cout << "Urban::ComputeTrueLength: tPath    773   //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength 
770   //         << " step= " << geomStepLength <<    774   //         << " step= " << geomStepLength << " &&& " << G4endl;
771                                                   775 
772   return tPathLength;                             776   return tPathLength;
773 }                                                 777 }
774                                                   778 
775 //....oooOO0OOooo........oooOO0OOooo........oo    779 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
776                                                   780 
777 G4ThreeVector&                                    781 G4ThreeVector& 
778 G4UrbanMscModel::SampleScattering(const G4Thre    782 G4UrbanMscModel::SampleScattering(const G4ThreeVector& oldDirection,
779                                   G4double /*s    783                                   G4double /*safety*/)
780 {                                                 784 {
781   fDisplacement.set(0.0,0.0,0.0);                 785   fDisplacement.set(0.0,0.0,0.0);
782   if(tPathLength >= currentRange) { return fDi    786   if(tPathLength >= currentRange) { return fDisplacement; }
783                                                   787 
784   G4double kinEnergy = currentKinEnergy;          788   G4double kinEnergy = currentKinEnergy;
785   if (tPathLength > currentRange*dtrl) {          789   if (tPathLength > currentRange*dtrl) {
786     kinEnergy = GetEnergy(particle,currentRang    790     kinEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
787   } else if(tPathLength > currentRange*0.01) {    791   } else if(tPathLength > currentRange*0.01) {
788     kinEnergy -= tPathLength*GetDEDX(particle,    792     kinEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple,
789                                      currentLo    793                                      currentLogKinEnergy);
790   }                                               794   }
791                                                   795 
792   if((tPathLength <= tlimitminfix) || (tPathLe    796   if((tPathLength <= tlimitminfix) || (tPathLength < tausmall*lambda0) || 
793      (kinEnergy <= CLHEP::eV)) { return fDispl    797      (kinEnergy <= CLHEP::eV)) { return fDisplacement; }
794                                                   798 
795   G4double cth = SampleCosineTheta(tPathLength    799   G4double cth = SampleCosineTheta(tPathLength,kinEnergy);
796                                                   800 
797   // protection against 'bad' cth values          801   // protection against 'bad' cth values
798   if(std::abs(cth) >= 1.0) { return fDisplacem    802   if(std::abs(cth) >= 1.0) { return fDisplacement; } 
799                                                   803 
800   G4double sth = std::sqrt((1.0 - cth)*(1.0 +     804   G4double sth = std::sqrt((1.0 - cth)*(1.0 + cth));
801   G4double phi = CLHEP::twopi*rndmEngineMod->f    805   G4double phi = CLHEP::twopi*rndmEngineMod->flat();
802   G4ThreeVector newDirection(sth*std::cos(phi)    806   G4ThreeVector newDirection(sth*std::cos(phi),sth*std::sin(phi),cth);
803   newDirection.rotateUz(oldDirection);            807   newDirection.rotateUz(oldDirection);
804                                                   808 
805   fParticleChange->ProposeMomentumDirection(ne    809   fParticleChange->ProposeMomentumDirection(newDirection);
806   /*                                              810   /*
807   G4cout << "G4UrbanMscModel::SampleSecondarie    811   G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
808          << " sinTheta= " << sth << " safety(m    812          << " sinTheta= " << sth << " safety(mm)= " << safety
809          << " trueStep(mm)= " << tPathLength      813          << " trueStep(mm)= " << tPathLength
810          << " geomStep(mm)= " << zPathLength      814          << " geomStep(mm)= " << zPathLength
811          << G4endl;                               815          << G4endl;
812   */                                              816   */
813                                                   817 
814   if (latDisplasment && currentTau >= tausmall    818   if (latDisplasment && currentTau >= tausmall) {
815     if(dispAlg96) { SampleDisplacement(sth, ph    819     if(dispAlg96) { SampleDisplacement(sth, phi); }
816     else          { SampleDisplacementNew(cth,    820     else          { SampleDisplacementNew(cth, phi); }
817     fDisplacement.rotateUz(oldDirection);         821     fDisplacement.rotateUz(oldDirection);
818   }                                               822   }
819   return fDisplacement;                           823   return fDisplacement;
820 }                                                 824 }
821                                                   825 
822 //....oooOO0OOooo........oooOO0OOooo........oo    826 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
823                                                   827 
824 G4double G4UrbanMscModel::SampleCosineTheta(G4    828 G4double G4UrbanMscModel::SampleCosineTheta(G4double trueStepLength,
825                                             G4    829                                             G4double kinEnergy)
826 {                                                 830 {
827   G4double cth = 1.0;                             831   G4double cth = 1.0;
828   G4double tau = trueStepLength/lambda0;          832   G4double tau = trueStepLength/lambda0;
829                                                   833 
830   // mean tau value                               834   // mean tau value
831   if(currentKinEnergy != kinEnergy) {             835   if(currentKinEnergy != kinEnergy) {
832     G4double lambda1 = GetTransportMeanFreePat    836     G4double lambda1 = GetTransportMeanFreePath(particle, kinEnergy);
833     if(std::abs(lambda1 - lambda0) > lambda0*0    837     if(std::abs(lambda1 - lambda0) > lambda0*0.01 && lambda1 > 0.) {
834       tau = trueStepLength*G4Log(lambda0/lambd    838       tau = trueStepLength*G4Log(lambda0/lambda1)/(lambda0-lambda1);
835     }                                             839     }
836   }                                               840   }
837                                                   841 
838   currentTau = tau;                               842   currentTau = tau;
839   lambdaeff = trueStepLength/currentTau;          843   lambdaeff = trueStepLength/currentTau;
840   currentRadLength = couple->GetMaterial()->Ge    844   currentRadLength = couple->GetMaterial()->GetRadlen();
841                                                   845 
842   if (tau >= taubig) { cth = -1.+2.*rndmEngine    846   if (tau >= taubig) { cth = -1.+2.*rndmEngineMod->flat(); }
843   else if (tau >= tausmall) {                     847   else if (tau >= tausmall) {
844     static const G4double numlim = 0.01;          848     static const G4double numlim = 0.01;
845     static const G4double onethird = 1./3.;       849     static const G4double onethird = 1./3.;
846     if(tau < numlim) {                            850     if(tau < numlim) {
847       xmeanth = 1.0 - tau*(1.0 - 0.5*tau);        851       xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
848       x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*one    852       x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*onethird;
849     } else {                                      853     } else {
850       xmeanth = G4Exp(-tau);                      854       xmeanth = G4Exp(-tau);
851       x2meanth = (1.+2.*G4Exp(-2.5*tau))*oneth    855       x2meanth = (1.+2.*G4Exp(-2.5*tau))*onethird;
852     }                                             856     }
853                                                   857 
854     // too large step of low-energy particle      858     // too large step of low-energy particle
855     G4double relloss = 1. - kinEnergy/currentK    859     G4double relloss = 1. - kinEnergy/currentKinEnergy;
856     static const G4double rellossmax= 0.50;       860     static const G4double rellossmax= 0.50;
857     if(relloss > rellossmax) {                    861     if(relloss > rellossmax) {
858       return SimpleScattering();                  862       return SimpleScattering();
859     }                                             863     }
860     // is step extreme small ?                    864     // is step extreme small ?
861     G4bool extremesmallstep = false;              865     G4bool extremesmallstep = false;
862     G4double tsmall = std::min(tlimitmin,lambd    866     G4double tsmall = std::min(tlimitmin,lambdalimit);
863                                                   867 
864     G4double theta0;                              868     G4double theta0;
865     if(trueStepLength > tsmall) {                 869     if(trueStepLength > tsmall) {
866       theta0 = ComputeTheta0(trueStepLength,ki    870       theta0 = ComputeTheta0(trueStepLength,kinEnergy);
867     } else {                                      871     } else {
868       theta0 = std::sqrt(trueStepLength/tsmall    872       theta0 = std::sqrt(trueStepLength/tsmall)
869   *ComputeTheta0(tsmall,kinEnergy);               873   *ComputeTheta0(tsmall,kinEnergy);
870       extremesmallstep = true;                    874       extremesmallstep = true;
871     }                                             875     }
872                                                   876 
873     static const G4double onesixth = 1./6.;       877     static const G4double onesixth = 1./6.;
874     static const G4double one12th = 1./12.;       878     static const G4double one12th = 1./12.;
875     static const G4double theta0max = CLHEP::p    879     static const G4double theta0max = CLHEP::pi*onesixth;
876                                                   880 
877     // protection for very small angles           881     // protection for very small angles
878     G4double theta2 = theta0*theta0;              882     G4double theta2 = theta0*theta0;
879                                                   883 
880     if(theta2 < tausmall) { return cth; }         884     if(theta2 < tausmall) { return cth; }
881     if(theta0 > theta0max) { return SimpleScat    885     if(theta0 > theta0max) { return SimpleScattering(); }
882                                                   886 
883     G4double x = theta2*(1.0 - theta2*one12th)    887     G4double x = theta2*(1.0 - theta2*one12th);
884     if(theta2 > numlim) {                         888     if(theta2 > numlim) {
885       G4double sth = 2*std::sin(0.5*theta0);      889       G4double sth = 2*std::sin(0.5*theta0);
886       x = sth*sth;                                890       x = sth*sth;
887     }                                             891     }
888                                                   892 
889     // parameter for tail                         893     // parameter for tail
890     G4double ltau = G4Log(tau);                   894     G4double ltau = G4Log(tau);
891     G4double u = !extremesmallstep ? G4Exp(lta    895     G4double u = !extremesmallstep ? G4Exp(ltau*onesixth) 
892       : G4Exp(G4Log(tsmall/lambda0)*onesixth);    896       : G4Exp(G4Log(tsmall/lambda0)*onesixth); 
893                                                   897 
894     G4double xx  = G4Log(lambdaeff/currentRadL    898     G4double xx  = G4Log(lambdaeff/currentRadLength);
895     G4double xsi = msc[idx]->coeffc1 +            899     G4double xsi = msc[idx]->coeffc1 + 
896       u*(msc[idx]->coeffc2+msc[idx]->coeffc3*u    900       u*(msc[idx]->coeffc2+msc[idx]->coeffc3*u)+msc[idx]->coeffc4*xx;
897                                                   901 
898     // tail should not be too big                 902     // tail should not be too big
899     xsi = std::max(xsi, 1.9);                     903     xsi = std::max(xsi, 1.9); 
900       /*                                          904       /*
901       if(KineticEnergy > 20*MeV && xsi < 1.6)     905       if(KineticEnergy > 20*MeV && xsi < 1.6) {
902         G4cout << "G4UrbanMscModel::SampleCosi    906         G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= " 
903                << KineticEnergy/GeV               907                << KineticEnergy/GeV 
904                << " !!** c= " << xsi              908                << " !!** c= " << xsi
905                << " **!! length(mm)= " << true    909                << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff 
906                << " " << couple->GetMaterial()    910                << " " << couple->GetMaterial()->GetName()
907                << " tau= " << tau << G4endl;      911                << " tau= " << tau << G4endl;
908       }                                           912       }
909       */                                          913       */
910                                                   914 
911     G4double c = xsi;                             915     G4double c = xsi;
912                                                   916 
913     if(std::abs(c-3.) < 0.001)      { c = 3.00    917     if(std::abs(c-3.) < 0.001)      { c = 3.001; }
914     else if(std::abs(c-2.) < 0.001) { c = 2.00    918     else if(std::abs(c-2.) < 0.001) { c = 2.001; }
915                                                   919 
916     G4double c1 = c-1.;                           920     G4double c1 = c-1.;
917     G4double ea = G4Exp(-xsi);                    921     G4double ea = G4Exp(-xsi);
918     G4double eaa = 1.-ea ;                        922     G4double eaa = 1.-ea ;
919     G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/ea    923     G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
920     G4double x0 = 1. - xsi*x;                     924     G4double x0 = 1. - xsi*x;
921                                                   925 
922     // G4cout << " xmean1= " << xmean1 << "  x    926     // G4cout << " xmean1= " << xmean1 << "  xmeanth= " << xmeanth << G4endl;
923                                                   927 
924     if(xmean1 <= 0.999*xmeanth) { return Simpl    928     if(xmean1 <= 0.999*xmeanth) { return SimpleScattering(); }
925                                                   929 
926     //from continuity of derivatives              930     //from continuity of derivatives
927     G4double b = 1.+(c-xsi)*x;                    931     G4double b = 1.+(c-xsi)*x;
928                                                   932 
929     G4double b1 = b+1.;                           933     G4double b1 = b+1.;
930     G4double bx = c*x;                            934     G4double bx = c*x;
931                                                   935 
932     G4double eb1 = G4Exp(G4Log(b1)*c1);           936     G4double eb1 = G4Exp(G4Log(b1)*c1);
933     G4double ebx = G4Exp(G4Log(bx)*c1);           937     G4double ebx = G4Exp(G4Log(bx)*c1);
934     G4double d = ebx/eb1;                         938     G4double d = ebx/eb1;
935                                                   939 
936     G4double xmean2 = (x0 + d - (bx - b1*d)/(c    940     G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
937                                                   941 
938     G4double f1x0 = ea/eaa;                       942     G4double f1x0 = ea/eaa;
939     G4double f2x0 = c1/(c*(1. - d));              943     G4double f2x0 = c1/(c*(1. - d));
940     G4double prob = f2x0/(f1x0+f2x0);             944     G4double prob = f2x0/(f1x0+f2x0);
941                                                   945 
942     G4double qprob = xmeanth/(prob*xmean1+(1.-    946     G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
943                                                   947 
944     // sampling of costheta                       948     // sampling of costheta
945     //G4cout << "c= " << c << " qprob= " << qp    949     //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
946     // << " c1= " << c1 << " b1= " << b1 << "     950     // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
947     //             << G4endl;                     951     //             << G4endl;
948     rndmEngineMod->flatArray(2, rndmarray);       952     rndmEngineMod->flatArray(2, rndmarray);
949     if(rndmarray[0] < qprob)                      953     if(rndmarray[0] < qprob)
950     {                                             954     {
951       G4double var = 0;                           955       G4double var = 0;
952       if(rndmarray[1] < prob) {                   956       if(rndmarray[1] < prob) {
953         cth = 1.+G4Log(ea+rndmEngineMod->flat(    957         cth = 1.+G4Log(ea+rndmEngineMod->flat()*eaa)*x;
954       } else {                                    958       } else {
955         var = (1.0 - d)*rndmEngineMod->flat();    959         var = (1.0 - d)*rndmEngineMod->flat();
956         if(var < numlim*d) {                      960         if(var < numlim*d) {
957           var /= (d*c1);                          961           var /= (d*c1); 
958           cth = -1.0 + var*(1.0 - 0.5*var*c)*(    962           cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
959         } else {                                  963         } else {
960           cth = 1. + x*(c - xsi - c*G4Exp(-G4L    964           cth = 1. + x*(c - xsi - c*G4Exp(-G4Log(var + d)/c1));
961         }                                         965         }
962       }                                           966       } 
963     } else {                                      967     } else {
964       cth = -1.+2.*rndmarray[1];                  968       cth = -1.+2.*rndmarray[1];
965     }                                             969     }
966   }                                               970   }
967   return cth;                                     971   return cth;
968 }                                                 972 }
969                                                   973 
970 //....oooOO0OOooo........oooOO0OOooo........oo    974 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
971                                                   975 
972 G4double G4UrbanMscModel::ComputeTheta0(G4doub    976 G4double G4UrbanMscModel::ComputeTheta0(G4double trueStepLength,
973                                         G4doub    977                                         G4double kinEnergy)
974 {                                                 978 {
975   // for all particles take the width of the c    979   // for all particles take the width of the central part
976   //  from a  parametrization similar to the H    980   //  from a  parametrization similar to the Highland formula
977   // ( Highland formula: Particle Physics Book    981   // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
978   G4double invbetacp = (kinEnergy+mass)/(kinEn    982   G4double invbetacp = (kinEnergy+mass)/(kinEnergy*(kinEnergy+2.*mass));
979   if(currentKinEnergy != kinEnergy) {             983   if(currentKinEnergy != kinEnergy) {
980     invbetacp = std::sqrt(invbetacp*(currentKi    984     invbetacp = std::sqrt(invbetacp*(currentKinEnergy+mass)/
981         (currentKinEnergy*(currentKinEnergy+2.    985         (currentKinEnergy*(currentKinEnergy+2.*mass)));
982   }                                               986   }
983   G4double y = trueStepLength/currentRadLength    987   G4double y = trueStepLength/currentRadLength;
984                                                   988 
985   if(fPosiCorrection && particle == positron)     989   if(fPosiCorrection && particle == positron)
986   {                                               990   {
987     static const G4double xl= 0.6;                991     static const G4double xl= 0.6;
988     static const G4double xh= 0.9;                992     static const G4double xh= 0.9;
989     static const G4double e = 113.0;              993     static const G4double e = 113.0;
990     G4double corr;                                994     G4double corr;
991                                                   995 
992     G4double tau = std::sqrt(currentKinEnergy*    996     G4double tau = std::sqrt(currentKinEnergy*kinEnergy)/mass;
993     G4double x = std::sqrt(tau*(tau+2.)/((tau+    997     G4double x = std::sqrt(tau*(tau+2.)/((tau+1.)*(tau+1.)));
994     G4double a = msc[idx]->posa;                  998     G4double a = msc[idx]->posa; 
995     G4double b = msc[idx]->posb;                  999     G4double b = msc[idx]->posb;
996     G4double c = msc[idx]->posc;                  1000     G4double c = msc[idx]->posc; 
997     G4double d = msc[idx]->posd;                  1001     G4double d = msc[idx]->posd;
998     if(x < xl) {                                  1002     if(x < xl) {
999       corr = a*(1.-G4Exp(-b*x));                  1003       corr = a*(1.-G4Exp(-b*x));  
1000     } else if(x > xh) {                          1004     } else if(x > xh) {
1001       corr = c+d*G4Exp(e*(x-1.));                1005       corr = c+d*G4Exp(e*(x-1.)); 
1002     } else {                                     1006     } else {
1003       G4double yl = a*(1.-G4Exp(-b*xl));         1007       G4double yl = a*(1.-G4Exp(-b*xl));
1004       G4double yh = c+d*G4Exp(e*(xh-1.));        1008       G4double yh = c+d*G4Exp(e*(xh-1.));
1005       G4double y0 = (yh-yl)/(xh-xl);             1009       G4double y0 = (yh-yl)/(xh-xl);
1006       G4double y1 = yl-y0*xl;                    1010       G4double y1 = yl-y0*xl;
1007       corr = y0*x+y1;                            1011       corr = y0*x+y1;
1008     }                                            1012     }
1009     //=======================================    1013     //==================================================================
1010     y *= corr*msc[idx]->pose;                    1014     y *= corr*msc[idx]->pose;
1011   }                                              1015   }
1012                                                  1016 
1013   static const G4double c_highland = 13.6*CLH    1017   static const G4double c_highland = 13.6*CLHEP::MeV;
1014   G4double theta0 = c_highland*std::abs(charg    1018   G4double theta0 = c_highland*std::abs(charge)*std::sqrt(y)*invbetacp;
1015                                                  1019  
1016   // correction factor from e- scattering dat    1020   // correction factor from e- scattering data
1017   theta0 *= (msc[idx]->coeffth1+msc[idx]->coe    1021   theta0 *= (msc[idx]->coeffth1+msc[idx]->coeffth2*G4Log(y));
1018   return theta0;                                 1022   return theta0;
1019 }                                                1023 }
1020                                                  1024 
1021 //....oooOO0OOooo........oooOO0OOooo........o    1025 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1022                                                  1026 
1023 void G4UrbanMscModel::SampleDisplacement(G4do    1027 void G4UrbanMscModel::SampleDisplacement(G4double, G4double phi)
1024 {                                                1028 { 
1025   // simple and fast sampling                    1029   // simple and fast sampling
1026   // based on single scattering results          1030   // based on single scattering results
1027   // u = r/rmax : mean value                     1031   // u = r/rmax : mean value
1028                                                  1032 
1029   G4double rmax = std::sqrt((tPathLength-zPat    1033   G4double rmax = std::sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1030   if(rmax > 0.)                                  1034   if(rmax > 0.)
1031   {                                              1035   {
1032     G4double r = 0.73*rmax;                      1036     G4double r = 0.73*rmax;
1033                                                  1037 
1034     // simple distribution for v=Phi-phi=psi     1038     // simple distribution for v=Phi-phi=psi ~exp(-beta*v)
1035     // beta determined from the requirement t    1039     // beta determined from the requirement that distribution should give
1036     // the same mean value than that obtained    1040     // the same mean value than that obtained from the ss simulation
1037                                                  1041 
1038     static const G4double cbeta  = 2.160;        1042     static const G4double cbeta  = 2.160;
1039     static const G4double cbeta1 = 1. - G4Exp    1043     static const G4double cbeta1 = 1. - G4Exp(-cbeta*CLHEP::pi);
1040     rndmEngineMod->flatArray(2, rndmarray);      1044     rndmEngineMod->flatArray(2, rndmarray);
1041     G4double psi = -G4Log(1. - rndmarray[0]*c    1045     G4double psi = -G4Log(1. - rndmarray[0]*cbeta1)/cbeta;
1042     G4double Phi = (rndmarray[1] < 0.5) ? phi    1046     G4double Phi = (rndmarray[1] < 0.5) ? phi+psi : phi-psi;
1043     fDisplacement.set(r*std::cos(Phi),r*std::    1047     fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1044   }                                              1048   }
1045 }                                                1049 }
1046                                                  1050 
1047 //....oooOO0OOooo........oooOO0OOooo........o    1051 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1048                                                  1052    
1049 void G4UrbanMscModel::SampleDisplacementNew(G    1053 void G4UrbanMscModel::SampleDisplacementNew(G4double, G4double phi)
1050 {                                                1054 {
1051   // simple and fast sampling                    1055   // simple and fast sampling
1052   // based on single scattering results          1056   // based on single scattering results
1053   // u = (r/rmax)**2 : distribution from ss s    1057   // u = (r/rmax)**2 : distribution from ss simulations
1054   const G4double eps = 1.e-3;                    1058   const G4double eps = 1.e-3;
1055   const G4double rmax =                          1059   const G4double rmax =
1056     std::sqrt((tPathLength-zPathLength)*(tPat    1060     std::sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1057                                                  1061   
1058   if(rmax > 0.)                                  1062   if(rmax > 0.)
1059   {                                              1063   {
1060     const G4double x0 = 0.73 ;                   1064     const G4double x0 = 0.73 ;
1061     const G4double alpha = G4Log(7.33)/x0 ;      1065     const G4double alpha = G4Log(7.33)/x0 ;
1062     const G4double a1 = 1.-x0 ;                  1066     const G4double a1 = 1.-x0 ;
1063     const G4double a2 = 1.-G4Exp(-alpha*x0) ;    1067     const G4double a2 = 1.-G4Exp(-alpha*x0) ;
1064     const G4double a3 = G4Exp(alpha*x0)-1. ;     1068     const G4double a3 = G4Exp(alpha*x0)-1. ;
1065     const G4double w1 = 2.*a2/(alpha*a1+2.*a2    1069     const G4double w1 = 2.*a2/(alpha*a1+2.*a2) ;
1066                                                  1070 
1067     G4double r, sqx;                             1071     G4double r, sqx;
1068     if (rmax/currentRange < eps)                 1072     if (rmax/currentRange < eps)
1069     {                                            1073     {
1070       r = 0.73*rmax ;                            1074       r = 0.73*rmax ;
1071       sqx = 1.;                                  1075       sqx = 1.;
1072     }                                            1076     }
1073     else                                         1077     else
1074     {                                            1078     {
1075       rndmEngineMod->flatArray(2,rndmarray);     1079       rndmEngineMod->flatArray(2,rndmarray);
1076       const G4double x = (rndmarray[0] < w1)     1080       const G4double x = (rndmarray[0] < w1) ? G4Log(1. + a3*rndmarray[1])/alpha :
1077         1. - a1*std::sqrt(1.-rndmarray[1]);      1081         1. - a1*std::sqrt(1.-rndmarray[1]);
1078                                                  1082 
1079       sqx = std::sqrt(x);                        1083       sqx = std::sqrt(x);
1080       r = sqx*rmax;                              1084       r = sqx*rmax;  
1081     }                                            1085     }
1082     // Gaussian distribution for Phi-phi=psi     1086     // Gaussian distribution for Phi-phi=psi 
1083     const G4double sigma = 0.1+0.9*sqx;          1087     const G4double sigma = 0.1+0.9*sqx;  
1084     const G4double psi = G4RandGauss::shoot(0    1088     const G4double psi = G4RandGauss::shoot(0.,sigma);
1085     const G4double Phi = phi+psi;                1089     const G4double Phi = phi+psi;  
1086     fDisplacement.set(r*std::cos(Phi), r*std:    1090     fDisplacement.set(r*std::cos(Phi), r*std::sin(Phi), 0.0);
1087   }                                              1091   }
1088 }                                                1092 }
1089                                                  1093 
1090 //....oooOO0OOooo........oooOO0OOooo........o    1094 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1091                                                  1095 
1092 void G4UrbanMscModel::InitialiseModelCache()     1096 void G4UrbanMscModel::InitialiseModelCache()                                   
1093 {                                                1097 {
1094   // it is assumed, that for the second run o    1098   // it is assumed, that for the second run only addition
1095   // of a new G4MaterialCutsCouple is possibl    1099   // of a new G4MaterialCutsCouple is possible
1096   auto theCoupleTable = G4ProductionCutsTable    1100   auto theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
1097   std::size_t numOfCouples = theCoupleTable->    1101   std::size_t numOfCouples = theCoupleTable->GetTableSize();
1098   if(numOfCouples != msc.size()) { msc.resize    1102   if(numOfCouples != msc.size()) { msc.resize(numOfCouples, nullptr); }
1099                                                  1103   
1100   for(G4int j=0; j<(G4int)numOfCouples; ++j)     1104   for(G4int j=0; j<(G4int)numOfCouples; ++j) {
1101     auto aCouple = theCoupleTable->GetMateria    1105     auto aCouple = theCoupleTable->GetMaterialCutsCouple(j);
1102                                                  1106 
1103     // new couple                                1107     // new couple
1104     msc[j] = new mscData();                      1108     msc[j] = new mscData();
1105     G4double Zeff = aCouple->GetMaterial()->G    1109     G4double Zeff = aCouple->GetMaterial()->GetIonisation()->GetZeffective();
1106     G4double sqrz = std::sqrt(Zeff);             1110     G4double sqrz = std::sqrt(Zeff);
1107     msc[j]->sqrtZ = sqrz;                        1111     msc[j]->sqrtZ = sqrz;
1108     // parameterisation of step limitation       1112     // parameterisation of step limitation
1109     msc[j]->factmin = dispAlg96 ? 0.001 : 0.0    1113     msc[j]->factmin = dispAlg96 ? 0.001 : 0.001/(1.+0.028*sqrz); 
1110     G4double lnZ = G4Log(Zeff);                  1114     G4double lnZ = G4Log(Zeff);
1111     // correction in theta0 formula              1115     // correction in theta0 formula
1112     G4double w = G4Exp(lnZ/6.);                  1116     G4double w = G4Exp(lnZ/6.);
1113     G4double facz = 0.990395+w*(-0.168386+w*0    1117     G4double facz = 0.990395+w*(-0.168386+w*0.093286);
1114     msc[j]->coeffth1 = facz*(1. - 8.7780e-2/Z    1118     msc[j]->coeffth1 = facz*(1. - 8.7780e-2/Zeff);
1115     msc[j]->coeffth2 = facz*(4.0780e-2 + 1.73    1119     msc[j]->coeffth2 = facz*(4.0780e-2 + 1.7315e-4*Zeff);
1116                                                  1120 
1117     // tail parameters                           1121     // tail parameters
1118     G4double Z13 = w*w;                          1122     G4double Z13 = w*w;
1119     msc[j]->coeffc1 = 2.3785    - Z13*(4.1981    1123     msc[j]->coeffc1 = 2.3785    - Z13*(4.1981e-1 - Z13*6.3100e-2);
1120     msc[j]->coeffc2 = 4.7526e-1 + Z13*(1.7694    1124     msc[j]->coeffc2 = 4.7526e-1 + Z13*(1.7694    - Z13*3.3885e-1);
1121     msc[j]->coeffc3 = 2.3683e-1 - Z13*(1.8111    1125     msc[j]->coeffc3 = 2.3683e-1 - Z13*(1.8111    - Z13*3.2774e-1);
1122     msc[j]->coeffc4 = 1.7888e-2 + Z13*(1.9659    1126     msc[j]->coeffc4 = 1.7888e-2 + Z13*(1.9659e-2 - Z13*2.6664e-3);
1123                                                  1127 
1124     msc[j]->Z23 = Z13*Z13;                       1128     msc[j]->Z23 = Z13*Z13;       
1125                                                  1129 
1126     msc[j]->stepmina = 27.725/(1.+0.203*Zeff)    1130     msc[j]->stepmina = 27.725/(1.+0.203*Zeff);
1127     msc[j]->stepminb =  6.152/(1.+0.111*Zeff)    1131     msc[j]->stepminb =  6.152/(1.+0.111*Zeff);
1128                                                  1132 
1129     // 21.07.2020                                1133     // 21.07.2020
1130     msc[j]->doverra = 9.6280e-1 - 8.4848e-2*m    1134     msc[j]->doverra = 9.6280e-1 - 8.4848e-2*msc[j]->sqrtZ + 4.3769e-3*Zeff;
1131                                                  1135 
1132     // 06.10.2020                                1136     // 06.10.2020 
1133     // msc[j]->doverra = 7.7024e-1 - 6.7878e-    1137     // msc[j]->doverra = 7.7024e-1 - 6.7878e-2*msc[j]->sqrtZ + 3.5015e-3*Zeff;
1134     msc[j]->doverrb = 1.15 - 9.76e-4*Zeff;       1138     msc[j]->doverrb = 1.15 - 9.76e-4*Zeff;
1135                                                  1139 
1136     // corrections for e+                        1140     // corrections for e+
1137     msc[j]->posa = 0.994-4.08e-3*Zeff;           1141     msc[j]->posa = 0.994-4.08e-3*Zeff;
1138     msc[j]->posb = 7.16+(52.6+365./Zeff)/Zeff    1142     msc[j]->posb = 7.16+(52.6+365./Zeff)/Zeff;
1139     msc[j]->posc = 1.000-4.47e-3*Zeff;           1143     msc[j]->posc = 1.000-4.47e-3*Zeff;
1140     msc[j]->posd = 1.21e-3*Zeff;                 1144     msc[j]->posd = 1.21e-3*Zeff;
1141     msc[j]->pose = 1.+Zeff*(1.84035e-4*Zeff-1    1145     msc[j]->pose = 1.+Zeff*(1.84035e-4*Zeff-1.86427e-2)+0.41125; 
1142   }                                              1146   }
1143 }                                                1147 }
1144                                                  1148 
1145 //....oooOO0OOooo........oooOO0OOooo........o    1149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1146                                                  1150