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 10.1.p3)


  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 // $Id:   $
                                                   >>  27 // GEANT4 tag $Name:  $
                                                   >>  28 //
 26 // -------------------------------------------     29 // -------------------------------------------------------------------
 27 //                                                 30 //   
 28 // GEANT4 Class file                               31 // GEANT4 Class file
 29 //                                                 32 //    
 30 //                                                 33 //
 31 // File name:   G4UrbanMscModel                    34 // File name:   G4UrbanMscModel
 32 //                                                 35 //
 33 // Author:      Laszlo Urban                       36 // Author:      Laszlo Urban
 34 //                                                 37 //
 35 // Creation date: 19.02.2013                       38 // Creation date: 19.02.2013
 36 //                                                 39 //
 37 // Created from G4UrbanMscModel96                  40 // Created from G4UrbanMscModel96
 38 //                                                 41 //
 39 // New parametrization for theta0                  42 // New parametrization for theta0
 40 // Correction for very small step length           43 // Correction for very small step length
 41 //                                                 44 //
 42 // Class Description:                              45 // Class Description:
 43 //                                                 46 //
 44 // Implementation of the model of multiple sca     47 // Implementation of the model of multiple scattering based on
 45 // H.W.Lewis Phys Rev 78 (1950) 526 and others     48 // H.W.Lewis Phys Rev 78 (1950) 526 and others
 46                                                    49 
 47 // -------------------------------------------     50 // -------------------------------------------------------------------
 48 // In its present form the model can be  used      51 // In its present form the model can be  used for simulation 
 49 //   of the e-/e+ multiple scattering              52 //   of the e-/e+ multiple scattering
                                                   >>  53 //
 50                                                    54 
 51                                                    55 
 52 //....oooOO0OOooo........oooOO0OOooo........oo     56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 53 //....oooOO0OOooo........oooOO0OOooo........oo     57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 54                                                    58 
 55 #include "G4UrbanMscModel.hh"                      59 #include "G4UrbanMscModel.hh"
 56 #include "G4PhysicalConstants.hh"                  60 #include "G4PhysicalConstants.hh"
 57 #include "G4SystemOfUnits.hh"                      61 #include "G4SystemOfUnits.hh"
 58 #include "Randomize.hh"                            62 #include "Randomize.hh"
 59 #include "G4Positron.hh"                       <<  63 #include "G4Electron.hh"
 60 #include "G4EmParameters.hh"                   <<  64 #include "G4LossTableManager.hh"
 61 #include "G4ParticleChangeForMSC.hh"               65 #include "G4ParticleChangeForMSC.hh"
 62 #include "G4ProductionCutsTable.hh"            << 
 63                                                    66 
 64 #include "G4Poisson.hh"                            67 #include "G4Poisson.hh"
 65 #include "G4Pow.hh"                                68 #include "G4Pow.hh"
                                                   >>  69 #include "globals.hh"
 66 #include "G4Log.hh"                                70 #include "G4Log.hh"
 67 #include "G4Exp.hh"                                71 #include "G4Exp.hh"
 68 #include "G4AutoLock.hh"                       << 
 69                                                    72 
 70 //....oooOO0OOooo........oooOO0OOooo........oo     73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 71                                                    74 
 72 std::vector<G4UrbanMscModel::mscData*> G4Urban <<  75 using namespace std;
 73                                                    76 
 74 namespace                                      <<  77 static const G4double Tlim = 10.*CLHEP::MeV;
 75 {                                              <<  78 static const G4double sigmafactor =
 76   G4Mutex theUrbanMutex = G4MUTEX_INITIALIZER; <<  79        CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
 77 }                                              <<  80 static const G4double epsfactor = 2.*CLHEP::electron_mass_c2*
                                                   >>  81        CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
                                                   >>  82        /(CLHEP::hbarc*CLHEP::hbarc);
                                                   >>  83 static const G4double beta2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
                                                   >>  84        ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHEP::electron_mass_c2));
                                                   >>  85 static const G4double bg2lim   = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
                                                   >>  86        (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
                                                   >>  87 
                                                   >>  88 static const G4double sig0[15] = {
                                                   >>  89  0.2672*CLHEP::barn,  0.5922*CLHEP::barn,  2.653*CLHEP::barn, 6.235*CLHEP::barn,
                                                   >>  90  11.69*CLHEP::barn  , 13.24*CLHEP::barn  , 16.12*CLHEP::barn, 23.00*CLHEP::barn,
                                                   >>  91  35.13*CLHEP::barn  , 39.95*CLHEP::barn  , 50.85*CLHEP::barn, 67.19*CLHEP::barn,
                                                   >>  92  91.15*CLHEP::barn  , 104.4*CLHEP::barn  , 113.1*CLHEP::barn};
                                                   >>  93 
                                                   >>  94 static const G4double Tdat[22] = { 
                                                   >>  95  100*CLHEP::eV,  200*CLHEP::eV,  400*CLHEP::eV,  700*CLHEP::eV,
                                                   >>  96    1*CLHEP::keV,   2*CLHEP::keV,   4*CLHEP::keV,   7*CLHEP::keV,
                                                   >>  97   10*CLHEP::keV,  20*CLHEP::keV,  40*CLHEP::keV,  70*CLHEP::keV,
                                                   >>  98  100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP::keV, 700*CLHEP::keV,
                                                   >>  99    1*CLHEP::MeV,   2*CLHEP::MeV,   4*CLHEP::MeV,   7*CLHEP::MeV,
                                                   >> 100   10*CLHEP::MeV,  20*CLHEP::MeV};
 78                                                   101 
 79 //....oooOO0OOooo........oooOO0OOooo........oo    102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 80                                                   103 
 81 G4UrbanMscModel::G4UrbanMscModel(const G4Strin    104 G4UrbanMscModel::G4UrbanMscModel(const G4String& nam)
 82   : G4VMscModel(nam)                              105   : G4VMscModel(nam)
 83 {                                                 106 {
 84   masslimite    = 0.6*CLHEP::MeV;              << 107   masslimite    = 0.6*MeV;
                                                   >> 108   lambdalimit   = 1.*mm;
 85   fr            = 0.02;                           109   fr            = 0.02;
 86   taubig        = 8.0;                            110   taubig        = 8.0;
 87   tausmall      = 1.e-16;                         111   tausmall      = 1.e-16;
 88   taulim        = 1.e-6;                          112   taulim        = 1.e-6;
 89   currentTau    = taulim;                         113   currentTau    = taulim;
 90   tlimitminfix  = 0.01*CLHEP::nm;              << 114   tlimitminfix  = 0.01*nm;             
 91   tlimitminfix2 = 1.*CLHEP::nm;                << 115   tlimitminfix2 =   1.*nm;             
 92   stepmin       = tlimitminfix;                   116   stepmin       = tlimitminfix;
 93   smallstep     = 1.e10;                          117   smallstep     = 1.e10;
 94   currentRange  = 0.;                          << 118   currentRange  = 0. ;
 95   rangeinit     = 0.;                             119   rangeinit     = 0.;
 96   tlimit        = 1.e10*CLHEP::mm;             << 120   tlimit        = 1.e10*mm;
 97   tlimitmin     = 10.*tlimitminfix;               121   tlimitmin     = 10.*tlimitminfix;            
 98   tgeom         = 1.e50*CLHEP::mm;             << 122   tgeom         = 1.e50*mm;
 99   geombig       = tgeom;                       << 123   geombig       = 1.e50*mm;
100   geommin       = 1.e-3*CLHEP::mm;             << 124   geommin       = 1.e-3*mm;
101   geomlimit     = geombig;                        125   geomlimit     = geombig;
102   presafety     = 0.;                          << 126   presafety     = 0.*mm;
103                                                   127 
104   positron      = G4Positron::Positron();      << 128   y             = 0.;
105   rndmEngineMod = G4Random::getTheEngine();    << 
106                                                   129 
107   drr = 0.35;                                  << 130   Zold          = 0.;
108   finalr = 10.*CLHEP::um;                      << 131   Zeff          = 1.;
109                                                << 132   Z2            = 1.;                
110   tlow = 5.*CLHEP::keV;                        << 133   Z23           = 1.;                    
111   invmev = 1.0/CLHEP::MeV;                     << 134   lnZ           = 0.;
                                                   >> 135   coeffth1      = 0.;
                                                   >> 136   coeffth2      = 0.;
                                                   >> 137   coeffc1       = 0.;
                                                   >> 138   coeffc2       = 0.;
                                                   >> 139   coeffc3       = 0.;
                                                   >> 140   coeffc4       = 0.;
                                                   >> 141 
                                                   >> 142   theta0max     = pi/6.;
                                                   >> 143   rellossmax    = 0.50;
                                                   >> 144   third         = 1./3.;
                                                   >> 145   particle      = 0;
                                                   >> 146   theManager    = G4LossTableManager::Instance(); 
                                                   >> 147   firstStep     = true; 
                                                   >> 148   inside        = false;  
                                                   >> 149   insideskin    = false;
                                                   >> 150   latDisplasmentbackup = false;
                                                   >> 151 
                                                   >> 152   rangecut = geombig;
                                                   >> 153   drr      = 0.35 ;
                                                   >> 154   finalr   = 10.*um ;
112                                                   155 
113   skindepth = skin*stepmin;                       156   skindepth = skin*stepmin;
114                                                   157 
115   mass = CLHEP::proton_mass_c2;                << 158   mass = proton_mass_c2;
116   charge = chargeSquare = 1.0;                 << 159   charge = ChargeSquare = 1.0;
117   currentKinEnergy = currentRadLength = lambda    160   currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength 
118     = zPathLength = par1 = par2 = par3 = rndma << 161     = zPathLength = par1 = par2 = par3 = 0;
119   currentLogKinEnergy = LOG_EKIN_MIN;          << 162 
                                                   >> 163   currentMaterialIndex = -1;
                                                   >> 164   fParticleChange = 0;
                                                   >> 165   couple = 0;
120 }                                                 166 }
121                                                   167 
122 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
123                                                   169 
124 G4UrbanMscModel::~G4UrbanMscModel()               170 G4UrbanMscModel::~G4UrbanMscModel()
125 {                                              << 171 {}
126   if(isFirstInstance) {                        << 
127     for(auto const & ptr : msc) { delete ptr;  << 
128     msc.clear();                               << 
129   }                                            << 
130 }                                              << 
131                                                   172 
132 //....oooOO0OOooo........oooOO0OOooo........oo    173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133                                                   174 
134 void G4UrbanMscModel::Initialise(const G4Parti    175 void G4UrbanMscModel::Initialise(const G4ParticleDefinition* p,
135                                  const G4DataV << 176          const G4DataVector&)
136 {                                                 177 {
                                                   >> 178   skindepth = skin*stepmin;
                                                   >> 179 
137   // set values of some data members              180   // set values of some data members
138   SetParticle(p);                                 181   SetParticle(p);
                                                   >> 182   /*
                                                   >> 183   if(p->GetPDGMass() > MeV) {
                                                   >> 184     G4cout << "### WARNING: G4UrbanMscModel model is used for " 
                                                   >> 185      << p->GetParticleName() << " !!! " << G4endl;
                                                   >> 186     G4cout << "###          This model should be used only for e+-" 
                                                   >> 187      << G4endl;
                                                   >> 188   }
                                                   >> 189   */
139   fParticleChange = GetParticleChangeForMSC(p)    190   fParticleChange = GetParticleChangeForMSC(p);
140   InitialiseParameters(p);                     << 
141                                                   191 
142   latDisplasmentbackup = latDisplasment;          192   latDisplasmentbackup = latDisplasment;
143                                                   193 
144   // if model is locked parameters should be d << 194   //G4cout << "### G4UrbanMscModel::Initialise done!" << G4endl;
145   if(!IsLocked()) {                            << 
146     dispAlg96 = G4EmParameters::Instance()->La << 
147     fPosiCorrection = G4EmParameters::Instance << 
148   }                                            << 
149                                                << 
150   // initialise cache only once                << 
151   if(0 == msc.size()) {                        << 
152     G4AutoLock l(&theUrbanMutex);              << 
153     if(0 == msc.size()) {                      << 
154       isFirstInstance = true;                  << 
155       msc.resize(1, nullptr);                  << 
156     }                                          << 
157     l.unlock();                                << 
158   }                                            << 
159   // initialise cache for each new run         << 
160   if(isFirstInstance) { InitialiseModelCache() << 
161                                                << 
162   /*                                           << 
163   G4cout << "### G4UrbanMscModel::Initialise d << 
164    << p->GetParticleName() << " type= " << ste << 
165   G4cout << "    RangeFact= " << facrange << " << 
166    << " SafetyFact= " << facsafety << " Lambda << 
167    << G4endl;                                  << 
168   */                                           << 
169 }                                                 195 }
170                                                   196 
171 //....oooOO0OOooo........oooOO0OOooo........oo    197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
172                                                   198 
173 G4double G4UrbanMscModel::ComputeCrossSectionP    199 G4double G4UrbanMscModel::ComputeCrossSectionPerAtom( 
174                              const G4ParticleD    200                              const G4ParticleDefinition* part,
175                                    G4double ki << 201                                    G4double KineticEnergy,
176                                    G4double at << 202                                    G4double AtomicNumber,G4double,
177                                    G4double, G << 203            G4double, G4double)
178 {                                                 204 {
179   static const G4double epsmin = 1.e-4 , epsma    205   static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
180                                                   206 
181   static const G4double Zdat[15] = { 4.,  6.,     207   static const G4double Zdat[15] = { 4.,  6., 13., 20., 26., 29., 32., 38.,47.,
182                                      50., 56., << 208              50., 56., 64., 74., 79., 82. };
183                                                   209 
184   // corr. factors for e-/e+ lambda for T <= T    210   // corr. factors for e-/e+ lambda for T <= Tlim
185   static const G4double celectron[15][22] =       211   static const G4double celectron[15][22] =
186           {{1.125,1.072,1.051,1.047,1.047,1.05    212           {{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    213             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    214             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    215            {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    216             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    217             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    218            {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    219             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    220             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    221            {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    222             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    223             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    224            {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    225             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    226             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    227            {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    228             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    229             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    230            {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    231             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    232             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    233            {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    234             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    235             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    236            {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    237             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    238             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    239            {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    240             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    241             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    242            {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    243             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    244             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    245            {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    246             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    247             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    248            {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    249             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    250             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    251            {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    252             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    253             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    254            {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    255             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    256             0.933,0.930,0.933,0.936,0.939,0.949            }};
231                                                << 257       
232   static const G4double cpositron[15][22] = {     258   static const G4double cpositron[15][22] = {
233            {2.589,2.044,1.658,1.446,1.347,1.21    259            {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    260             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    261             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    262            {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    263             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    264             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    265            {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    266             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    267             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    268            {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    269             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    270             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    271            {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    272             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    273             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    274            {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    275             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    276             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    277            {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    278             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    279             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    280            {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    281             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    282             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    283            {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    284             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    285             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    286            {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    287             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    288             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    289            {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    290             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    291             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    292            {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    293             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    294             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    295            {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    296             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    297             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    298            {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    299             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    300             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    301            {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    302             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    303             1.456,1.412,1.364,1.328,1.307,1.282            }};
278                                                   304 
279   //data/corrections for T > Tlim                 305   //data/corrections for T > Tlim  
280                                                << 306                      
281   static const G4double hecorr[15] = {            307   static const G4double hecorr[15] = {
282     120.70, 117.50, 105.00, 92.92, 79.23,  74.    308     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    309     57.39,  41.97,  36.14, 24.53, 10.21,  -7.855, -16.84,
284     -22.30};                                      310     -22.30};
285                                                   311 
286   G4double sigma;                                 312   G4double sigma;
287   SetParticle(part);                              313   SetParticle(part);
288                                                   314 
289   G4double Z23 = G4Pow::GetInstance()->Z23(G4l << 315   Z23 = G4Pow::GetInstance()->Z23(G4lrint(AtomicNumber));
290                                                   316 
291   // correction if particle .ne. e-/e+            317   // correction if particle .ne. e-/e+
292   // compute equivalent kinetic energy            318   // compute equivalent kinetic energy
293   // lambda depends on p*beta ....                319   // lambda depends on p*beta ....
294                                                   320 
295   G4double eKineticEnergy = kinEnergy;         << 321   G4double eKineticEnergy = KineticEnergy;
296                                                   322 
297   if(mass > CLHEP::electron_mass_c2)           << 323   if(mass > electron_mass_c2)
298   {                                               324   {
299      G4double TAU = kinEnergy/mass ;           << 325      G4double TAU = KineticEnergy/mass ;
300      G4double c = mass*TAU*(TAU+2.)/(CLHEP::el << 326      G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
301      G4double w = c-2.;                        << 327      G4double w = c-2. ;
302      G4double tau = 0.5*(w+std::sqrt(w*w+4.*c) << 328      G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
303      eKineticEnergy = CLHEP::electron_mass_c2* << 329      eKineticEnergy = electron_mass_c2*tau ;
304   }                                               330   }
305                                                   331 
306   G4double eTotalEnergy = eKineticEnergy + CLH << 332   G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
307   G4double beta2 = eKineticEnergy*(eTotalEnerg << 333   G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
308                                  /(eTotalEnerg    334                                  /(eTotalEnergy*eTotalEnergy);
309   G4double bg2   = eKineticEnergy*(eTotalEnerg << 335   G4double bg2   = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
310                                  /(CLHEP::elec << 336                                  /(electron_mass_c2*electron_mass_c2);
311                                                   337 
312   static const G4double epsfactor = 2.*CLHEP:: << 
313     CLHEP::electron_mass_c2*CLHEP::Bohr_radius << 
314     /(CLHEP::hbarc*CLHEP::hbarc);              << 
315   G4double eps = epsfactor*bg2/Z23;               338   G4double eps = epsfactor*bg2/Z23;
316                                                   339 
317   if     (eps<epsmin)  sigma = 2.*eps*eps;        340   if     (eps<epsmin)  sigma = 2.*eps*eps;
318   else if(eps<epsmax)  sigma = G4Log(1.+2.*eps    341   else if(eps<epsmax)  sigma = G4Log(1.+2.*eps)-2.*eps/(1.+2.*eps);
319   else                 sigma = G4Log(2.*eps)-1    342   else                 sigma = G4Log(2.*eps)-1.+1./eps;
320                                                   343 
321   sigma *= chargeSquare*atomicNumber*atomicNum << 344   sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
322                                                   345 
323   // interpolate in AtomicNumber and beta2        346   // interpolate in AtomicNumber and beta2 
324   G4double c1,c2,cc1;                          << 347   G4double c1,c2,cc1,cc2,corr;
325                                                   348 
326   // get bin number in Z                          349   // get bin number in Z
327   G4int iZ = 14;                                  350   G4int iZ = 14;
328   // Loop checking, 03-Aug-2015, Vladimir Ivan << 351   while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
329   while ((iZ>=0)&&(Zdat[iZ]>=atomicNumber)) {  << 352   if (iZ==14)                               iZ = 13;
330                                                << 353   if (iZ==-1)                               iZ = 0 ;
331   iZ = std::min(std::max(iZ, 0), 13);          << 
332                                                   354 
333   G4double ZZ1 = Zdat[iZ];                        355   G4double ZZ1 = Zdat[iZ];
334   G4double ZZ2 = Zdat[iZ+1];                      356   G4double ZZ2 = Zdat[iZ+1];
335   G4double ratZ = (atomicNumber-ZZ1)*(atomicNu << 357   G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
336                   ((ZZ2-ZZ1)*(ZZ2+ZZ1));          358                   ((ZZ2-ZZ1)*(ZZ2+ZZ1));
337                                                   359 
338   static const G4double Tlim = 10.*CLHEP::MeV; << 
339   static const G4double sigmafactor =          << 
340     CLHEP::twopi*CLHEP::classic_electr_radius* << 
341   static const G4double beta2lim = Tlim*(Tlim+ << 
342     ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHE << 
343   static const G4double bg2lim   = Tlim*(Tlim+ << 
344     (CLHEP::electron_mass_c2*CLHEP::electron_m << 
345                                                << 
346   static const G4double sig0[15] = {           << 
347     0.2672*CLHEP::barn,  0.5922*CLHEP::barn,   << 
348     11.69*CLHEP::barn  , 13.24*CLHEP::barn  ,  << 
349     35.13*CLHEP::barn  , 39.95*CLHEP::barn  ,  << 
350     91.15*CLHEP::barn  , 104.4*CLHEP::barn  ,  << 
351                                                << 
352   static const G4double Tdat[22] = {           << 
353     100*CLHEP::eV,  200*CLHEP::eV,  400*CLHEP: << 
354     1*CLHEP::keV,   2*CLHEP::keV,   4*CLHEP::k << 
355     10*CLHEP::keV,  20*CLHEP::keV,  40*CLHEP:: << 
356     100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP: << 
357     1*CLHEP::MeV,   2*CLHEP::MeV,   4*CLHEP::M << 
358     10*CLHEP::MeV,  20*CLHEP::MeV};            << 
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 << 
365     while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)    364     while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
366                                                << 365     if(iT==21)                                  iT = 20;
367     iT = std::min(std::max(iT, 0), 20);        << 366     if(iT==-1)                                  iT = 0 ;
368                                                   367 
369     //  calculate betasquare values               368     //  calculate betasquare values
370     G4double T = Tdat[iT];                     << 369     G4double T = Tdat[iT],   E = T + electron_mass_c2;
371     G4double E = T + CLHEP::electron_mass_c2;  << 370     G4double b2small = T*(E+electron_mass_c2)/(E*E);
372     G4double b2small = T*(E+CLHEP::electron_ma << 371 
373                                                << 372     T = Tdat[iT+1]; E = T + electron_mass_c2;
374     T = Tdat[iT+1];                            << 373     G4double b2big = T*(E+electron_mass_c2)/(E*E);
375     E = T + CLHEP::electron_mass_c2;           << 
376     G4double b2big = T*(E+CLHEP::electron_mass << 
377     G4double ratb2 = (beta2-b2small)/(b2big-b2    374     G4double ratb2 = (beta2-b2small)/(b2big-b2small);
378                                                   375 
379     if (charge < 0.)                              376     if (charge < 0.)
380     {                                             377     {
381        c1 = celectron[iZ][iT];                    378        c1 = celectron[iZ][iT];
382        c2 = celectron[iZ+1][iT];                  379        c2 = celectron[iZ+1][iT];
383        cc1 = c1+ratZ*(c2-c1);                     380        cc1 = c1+ratZ*(c2-c1);
384                                                   381 
385        c1 = celectron[iZ][iT+1];                  382        c1 = celectron[iZ][iT+1];
386        c2 = celectron[iZ+1][iT+1];                383        c2 = celectron[iZ+1][iT+1];
                                                   >> 384        cc2 = c1+ratZ*(c2-c1);
                                                   >> 385 
                                                   >> 386        corr = cc1+ratb2*(cc2-cc1);
                                                   >> 387 
                                                   >> 388        sigma *= sigmafactor/corr;
387     }                                             389     }
388     else                                          390     else              
389     {                                             391     {
390        c1 = cpositron[iZ][iT];                    392        c1 = cpositron[iZ][iT];
391        c2 = cpositron[iZ+1][iT];                  393        c2 = cpositron[iZ+1][iT];
392        cc1 = c1+ratZ*(c2-c1);                     394        cc1 = c1+ratZ*(c2-c1);
393                                                   395 
394        c1 = cpositron[iZ][iT+1];                  396        c1 = cpositron[iZ][iT+1];
395        c2 = cpositron[iZ+1][iT+1];                397        c2 = cpositron[iZ+1][iT+1];
                                                   >> 398        cc2 = c1+ratZ*(c2-c1);
                                                   >> 399 
                                                   >> 400        corr = cc1+ratb2*(cc2-cc1);
                                                   >> 401 
                                                   >> 402        sigma *= sigmafactor/corr;
396     }                                             403     }
397     G4double cc2 = c1+ratZ*(c2-c1);            << 
398     sigma *= sigmafactor/(cc1+ratb2*(cc2-cc1)) << 
399   }                                               404   }
400   else                                            405   else
401   {                                               406   {
402     c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2    407     c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
403     c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(b    408     c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
404     if((atomicNumber >= ZZ1) && (atomicNumber  << 409     if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
405       sigma = c1+ratZ*(c2-c1) ;                   410       sigma = c1+ratZ*(c2-c1) ;
406     else if(atomicNumber < ZZ1)                << 411     else if(AtomicNumber < ZZ1)
407       sigma = atomicNumber*atomicNumber*c1/(ZZ << 412       sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
408     else if(atomicNumber > ZZ2)                << 413     else if(AtomicNumber > ZZ2)
409       sigma = atomicNumber*atomicNumber*c2/(ZZ << 414       sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
410   }                                               415   }
411   // low energy correction based on theory     << 
412   sigma *= (1.+0.30/(1.+std::sqrt(1000.*eKinet << 
413                                                << 
414   return sigma;                                   416   return sigma;
                                                   >> 417 
415 }                                                 418 }
416                                                   419 
417 //....oooOO0OOooo........oooOO0OOooo........oo    420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
418                                                   421 
419 void G4UrbanMscModel::StartTracking(G4Track* t    422 void G4UrbanMscModel::StartTracking(G4Track* track)
420 {                                                 423 {
421   SetParticle(track->GetDynamicParticle()->Get    424   SetParticle(track->GetDynamicParticle()->GetDefinition());
422   firstStep = true;                               425   firstStep = true; 
                                                   >> 426   inside = false;
423   insideskin = false;                             427   insideskin = false;
424   fr = facrange;                               << 428   tlimit = geombig;
425   tlimit = tgeom = rangeinit = geombig;        << 429   stepmin = tlimitminfix ;
426   smallstep     = 1.e10;                       << 430   tlimitmin = 10.*stepmin ;
427   stepmin       = tlimitminfix;                << 431   G4VEmModel::StartTracking(track);
428   tlimitmin     = 10.*tlimitminfix;            << 
429   rndmEngineMod = G4Random::getTheEngine();    << 
430 }                                                 432 }
431                                                   433 
432 //....oooOO0OOooo........oooOO0OOooo........oo    434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
433                                                   435 
434 G4double G4UrbanMscModel::ComputeTruePathLengt    436 G4double G4UrbanMscModel::ComputeTruePathLengthLimit(
435                              const G4Track& tr    437                              const G4Track& track,
436                              G4double& current << 438            G4double& currentMinimalStep)
437 {                                                 439 {
438   tPathLength = currentMinimalStep;               440   tPathLength = currentMinimalStep;
439   const G4DynamicParticle* dp = track.GetDynam    441   const G4DynamicParticle* dp = track.GetDynamicParticle();
440                                                   442   
441   G4StepPoint* sp = track.GetStep()->GetPreSte    443   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
442   G4StepStatus stepStatus = sp->GetStepStatus(    444   G4StepStatus stepStatus = sp->GetStepStatus();
443   couple = track.GetMaterialCutsCouple();         445   couple = track.GetMaterialCutsCouple();
444   SetCurrentCouple(couple);                       446   SetCurrentCouple(couple); 
445   idx = couple->GetIndex();                    << 447   currentMaterialIndex = couple->GetIndex();
446   currentKinEnergy = dp->GetKineticEnergy();      448   currentKinEnergy = dp->GetKineticEnergy();
447   currentLogKinEnergy = dp->GetLogKineticEnerg << 449   currentRange = GetRange(particle,currentKinEnergy,couple);
448   currentRange = GetRange(particle,currentKinE << 450   lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
449   lambda0 = GetTransportMeanFreePath(particle, << 451   tPathLength = min(tPathLength,currentRange);
450                                                << 452 
451   tPathLength = std::min(tPathLength,currentRa << 453   // set flag to default values
452   /*                                           << 454   latDisplasment = latDisplasmentbackup;
453   G4cout << "G4Urban::StepLimit tPathLength= " << 455   /*
454   << " range= " <<currentRange<< " lambda= "<< << 456   G4cout << "G4Urban::StepLimit tPathLength= " 
455             <<G4endl;                          << 457      <<tPathLength<<" inside= " << inside
                                                   >> 458    << " range= " <<currentRange<< " lambda= "<<lambda0
                                                   >> 459      <<G4endl;
456   */                                              460   */
457   // extreme small step                        << 461   // stop here if small range particle
458   if(tPathLength < tlimitminfix) {             << 462   if(inside) { 
459     latDisplasment = false;                    << 463     latDisplasment = false;   
460     return ConvertTrueToGeom(tPathLength, curr << 464     return ConvertTrueToGeom(tPathLength, currentMinimalStep); 
461   }                                               465   }
462                                                << 466   // stop here if small step
463   presafety = (stepStatus == fGeomBoundary) ?  << 467   if(tPathLength < tlimitminfix) { 
464               : ComputeSafety(sp->GetPosition( << 
465                                                << 
466   // stop here if small step or range is less  << 
467   if((tPathLength == currentRange && tPathLeng << 
468      tPathLength < tlimitminfix) {             << 
469     latDisplasment = false;                       468     latDisplasment = false;   
470     return ConvertTrueToGeom(tPathLength, curr    469     return ConvertTrueToGeom(tPathLength, currentMinimalStep); 
471   }                                               470   }
472                                                << 471   
473   // upper limit for the straight line distanc << 472   presafety = sp->GetSafety();
474   // for electrons and positrons               << 
475   G4double distance = (mass < masslimite)      << 
476     ? currentRange*msc[idx]->doverra           << 
477     // for muons, hadrons                      << 
478     : currentRange*msc[idx]->doverrb;          << 
479                                                << 
480   /*                                              473   /*  
481   G4cout << "G4Urban::StepLimit tPathLength= "    474   G4cout << "G4Urban::StepLimit tPathLength= " 
482             <<tPathLength<<" safety= " << pres << 475      <<tPathLength<<" safety= " << presafety
483           << " range= " <<currentRange<< " lam    476           << " range= " <<currentRange<< " lambda= "<<lambda0
484             << " Alg: " << steppingAlgorithm < << 477      << " Alg: " << steppingAlgorithm <<G4endl;
485   */                                              478   */
486   // far from geometry boundary                   479   // far from geometry boundary
487   if(distance < presafety)                     << 480   if(currentRange < presafety)
488     {                                             481     {
                                                   >> 482       inside = true;
489       latDisplasment = false;                     483       latDisplasment = false;   
490       return ConvertTrueToGeom(tPathLength, cu    484       return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
491     }                                             485     }
492                                                   486 
493   latDisplasment = latDisplasmentbackup;       << 487   // standard  version
494   // ----------------------------------------- << 488   //
495   // distance to boundary                      << 
496   if (steppingAlgorithm == fUseDistanceToBound    489   if (steppingAlgorithm == fUseDistanceToBoundary)
497     {                                             490     {
498       //compute geomlimit and presafety           491       //compute geomlimit and presafety 
499       geomlimit = ComputeGeomLimit(track, pres    492       geomlimit = ComputeGeomLimit(track, presafety, currentRange);
500       /*                                       << 493 
501         G4cout << "G4Urban::Distance to bounda << 494       // is it far from boundary ?
502             <<geomlimit<<" safety= " << presaf << 495       if(currentRange < presafety)
503       */                                       << 496   {
                                                   >> 497     inside = true;
                                                   >> 498           latDisplasment = false;   
                                                   >> 499     return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
                                                   >> 500   }
504                                                   501 
505       smallstep += 1.;                            502       smallstep += 1.;
506       insideskin = false;                         503       insideskin = false;
507       tgeom = geombig;                         << 
508                                                   504 
509       // initialisation at first step and at t << 
510       if(firstStep || (stepStatus == fGeomBoun    505       if(firstStep || (stepStatus == fGeomBoundary))
511         {                                         506         {
512           rangeinit = currentRange;               507           rangeinit = currentRange;
513           if(!firstStep) { smallstep = 1.; }   << 508           if(firstStep) { smallstep = 1.e10; }
                                                   >> 509           else          { smallstep = 1.; }
514                                                   510 
                                                   >> 511           //define stepmin here (it depends on lambda!)
                                                   >> 512           //rough estimation of lambda_elastic/lambda_transport
                                                   >> 513           G4double rat = currentKinEnergy/MeV ;
                                                   >> 514           rat = 1.e-3/(rat*(10.+rat)) ;
515           //stepmin ~ lambda_elastic              515           //stepmin ~ lambda_elastic
516           stepmin = ComputeStepmin();          << 516           stepmin = rat*lambda0;
517           skindepth = skin*stepmin;               517           skindepth = skin*stepmin;
518           tlimitmin = ComputeTlimitmin();      << 518           //define tlimitmin
519         /*                                     << 519           tlimitmin = 10.*stepmin;
520           G4cout << "rangeinit= " << rangeinit << 520           tlimitmin = max(tlimitmin,tlimitminfix);
521                  << " tlimitmin= " << tlimitmi << 521     //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
522                  << geomlimit <<G4endl;        << 522     //       << " tlimitmin= " << tlimitmin << " geomlimit= " 
523         */                                     << 523     //       << geomlimit <<G4endl;
524         }                                      << 524           // constraint from the geometry
525       // constraint from the geometry          << 525           if((geomlimit < geombig) && (geomlimit > geommin))
526       if((geomlimit < geombig) && (geomlimit > << 526             {
527         {                                      << 527               // geomlimit is a geometrical step length
528           // geomlimit is a geometrical step l << 528               // transform it to true path length (estimation)
529           // transform it to true path length  << 529               if((1.-geomlimit/lambda0) > 0.)
530           if(lambda0 > geomlimit) {            << 530                 geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin ;
531             geomlimit = -lambda0*G4Log(1.-geom << 531 
532           }                                    << 532               if(stepStatus == fGeomBoundary)
533           tgeom = (stepStatus == fGeomBoundary << 533                 tgeom = geomlimit/facgeom;
534       : facrange*rangeinit + stepmin;          << 534               else
                                                   >> 535                 tgeom = 2.*geomlimit/facgeom;
                                                   >> 536             }
                                                   >> 537             else
                                                   >> 538               tgeom = geombig;
535         }                                         539         }
536                                                   540 
537       //step limit                                541       //step limit 
538       tlimit = (currentRange > presafety) ?    << 542       tlimit = facrange*rangeinit;              
539         std::max(facrange*rangeinit, facsafety << 
540                                                   543 
541       //lower limit for tlimit                    544       //lower limit for tlimit
542       tlimit = std::min(std::max(tlimit,tlimit << 545       tlimit = max(tlimit,tlimitmin);
543       /*                                       << 546 
544       G4cout << "tgeom= " << tgeom << " geomli << 547       tlimit = min(tlimit,tgeom); 
545             << " tlimit= " << tlimit << " pres << 548 
546       */                                       << 549       //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit  
                                                   >> 550       //      << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
                                                   >> 551 
547       // shortcut                                 552       // shortcut
548       if((tPathLength < tlimit) && (tPathLengt    553       if((tPathLength < tlimit) && (tPathLength < presafety) &&
549          (smallstep > skin) && (tPathLength <     554          (smallstep > skin) && (tPathLength < geomlimit-0.999*skindepth))
550       {                                           555       {
551         return ConvertTrueToGeom(tPathLength,  << 556   return ConvertTrueToGeom(tPathLength, currentMinimalStep);   
552       }                                           557       }
553                                                   558 
554       // step reduction near to boundary          559       // step reduction near to boundary
555       if(smallstep <= skin)                       560       if(smallstep <= skin)
556         {                                      << 561   {
557           tlimit = stepmin;                    << 562     tlimit = stepmin;
558           insideskin = true;                   << 563     insideskin = true;
559         }                                      << 564   }
560       else if(geomlimit < geombig)                565       else if(geomlimit < geombig)
561         {                                      << 566   {
562           if(geomlimit > skindepth)            << 567     if(geomlimit > skindepth)
563             {                                  << 568       {
564               tlimit = std::min(tlimit, geomli << 569         tlimit = min(tlimit, geomlimit-0.999*skindepth);
565             }                                  << 570       }
566           else                                 << 571     else
567             {                                  << 572       {
568               insideskin = true;               << 573         insideskin = true;
569               tlimit = std::min(tlimit, stepmi << 574         tlimit = min(tlimit, stepmin);
570             }                                  << 575       }
571         }                                      << 576   }
                                                   >> 577 
                                                   >> 578       tlimit = max(tlimit, stepmin); 
                                                   >> 579 
                                                   >> 580       // randomize 1st step or 1st 'normal' step in volume
                                                   >> 581       if(firstStep || ((smallstep == skin+1) && !insideskin)) 
                                                   >> 582         { 
                                                   >> 583           G4double temptlimit = tlimit;
                                                   >> 584           if(temptlimit > tlimitmin)
                                                   >> 585           {
                                                   >> 586             do {
                                                   >> 587             temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.3*tlimit); 
                                                   >> 588       } while ((temptlimit < tlimitmin) || 
                                                   >> 589          (temptlimit > 2.*tlimit-tlimitmin));
                                                   >> 590           }
                                                   >> 591           else { temptlimit = tlimitmin; }
572                                                   592 
573       tlimit = std::max(tlimit, stepmin);      << 593           tPathLength = min(tPathLength, temptlimit);
                                                   >> 594         }
                                                   >> 595       else
                                                   >> 596         {  
                                                   >> 597           tPathLength = min(tPathLength, tlimit); 
                                                   >> 598         }
574                                                   599 
575       // randomise if not 'small' step and ste << 
576       tPathLength = (tlimit < tPathLength && s << 
577         ? std::min(tPathLength, Randomizetlimi << 
578   : std::min(tPathLength, tlimit);             << 
579     }                                             600     }
580   // ----------------------------------------- << 601     // for 'normal' simulation with or without magnetic field 
581   // for simulation with or without magnetic f << 602     //  there no small step/single scattering at boundaries
582   // there no small step/single scattering at  << 
583   else if(steppingAlgorithm == fUseSafety)        603   else if(steppingAlgorithm == fUseSafety)
584     {                                             604     {
585       if(firstStep || (stepStatus == fGeomBoun << 605       if(currentRange < presafety)
                                                   >> 606         {
                                                   >> 607           inside = true;
                                                   >> 608           latDisplasment = false;
                                                   >> 609           return ConvertTrueToGeom(tPathLength, currentMinimalStep);
                                                   >> 610         }   
                                                   >> 611       else if(stepStatus != fGeomBoundary)  {
                                                   >> 612   presafety = ComputeSafety(sp->GetPosition(),tPathLength); 
                                                   >> 613       }
                                                   >> 614       /*
                                                   >> 615       G4cout << "presafety= " << presafety
                                                   >> 616        << " firstStep= " << firstStep
                                                   >> 617        << " stepStatus= " << stepStatus 
                                                   >> 618        << G4endl;
                                                   >> 619       */
                                                   >> 620       // is far from boundary
                                                   >> 621       if(currentRange < presafety)
                                                   >> 622         {
                                                   >> 623           inside = true;
                                                   >> 624           latDisplasment = false;
                                                   >> 625           return ConvertTrueToGeom(tPathLength, currentMinimalStep);  
                                                   >> 626         }
                                                   >> 627 
                                                   >> 628       if(firstStep || stepStatus == fGeomBoundary)
                                                   >> 629       {
586         rangeinit = currentRange;                 630         rangeinit = currentRange;
587         fr = facrange;                            631         fr = facrange;
588         // stepping for e+/e- only (not for mu << 632         // 9.1 like stepping for e+/e- only (not for muons,hadrons)
589         if(mass < masslimite)                     633         if(mass < masslimite) 
590           {                                    << 634         {
591             rangeinit = std::max(rangeinit, la << 635     rangeinit = max(rangeinit, lambda0);
592             if(lambda0 > lambdalimit) {        << 636           if(lambda0 > lambdalimit) {
593               fr *= (0.75+0.25*lambda0/lambdal << 637             fr *= (0.75+0.25*lambda0/lambdalimit);
594             }                                  << 638     }
595           }                                    << 639         }
                                                   >> 640 
596         //lower limit for tlimit                  641         //lower limit for tlimit
597         stepmin = ComputeStepmin();            << 642         G4double rat = currentKinEnergy/MeV;
598         tlimitmin = ComputeTlimitmin();        << 643         rat = 1.e-3/(rat*(10 + rat)) ;
                                                   >> 644         stepmin = lambda0*rat;
                                                   >> 645         tlimitmin = max(10*stepmin, tlimitminfix);
599       }                                           646       }
600                                                   647 
601       //step limit                                648       //step limit
602       tlimit = (currentRange > presafety) ?    << 649       tlimit = max(fr*rangeinit, facsafety*presafety);
603         std::max(fr*rangeinit, facsafety*presa << 
604                                                   650   
605       //lower limit for tlimit                    651       //lower limit for tlimit
606       tlimit = std::max(tlimit, tlimitmin);    << 652       tlimit = max(tlimit, tlimitmin); 
607                                                << 653 
608       // randomise if step determined by msc   << 654       if(firstStep || stepStatus == fGeomBoundary)
609       tPathLength = (tlimit < tPathLength) ?   << 655       {
610         std::min(tPathLength, Randomizetlimit( << 656         G4double temptlimit = tlimit;
                                                   >> 657         if(temptlimit > tlimitmin)
                                                   >> 658         {
                                                   >> 659           do {
                                                   >> 660             temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.3*tlimit); 
                                                   >> 661     } while ((temptlimit < tlimitmin) ||
                                                   >> 662        (temptlimit > 2.*tlimit-tlimitmin));
                                                   >> 663         }
                                                   >> 664         else { temptlimit = tlimitmin; }
                                                   >> 665 
                                                   >> 666         tPathLength = min(tPathLength, temptlimit);
                                                   >> 667       }
                                                   >> 668       else { tPathLength = min(tPathLength, tlimit); }
611     }                                             669     }
612   // ----------------------------------------- << 670   // new stepping mode UseSafetyPlus
613   // for simulation with or without magnetic f << 671     else if(steppingAlgorithm == fUseSafetyPlus)
614   // there is small step/single scattering at  << 
615   else if(steppingAlgorithm == fUseSafetyPlus) << 
616     {                                             672     {
617       if(firstStep || (stepStatus == fGeomBoun << 673       if(currentRange < presafety)
                                                   >> 674         {
                                                   >> 675           inside = true;
                                                   >> 676           latDisplasment = false;
                                                   >> 677           return ConvertTrueToGeom(tPathLength, currentMinimalStep);
                                                   >> 678         }
                                                   >> 679       else if(stepStatus != fGeomBoundary)  {
                                                   >> 680         presafety = ComputeSafety(sp->GetPosition(),tPathLength);
                                                   >> 681       }
                                                   >> 682       /*
                                                   >> 683       G4cout << "presafety= " << presafety
                                                   >> 684              << " firstStep= " << firstStep
                                                   >> 685              << " stepStatus= " << stepStatus
                                                   >> 686              << G4endl;
                                                   >> 687       */
                                                   >> 688       // is far from boundary
                                                   >> 689       if(currentRange < presafety)
                                                   >> 690         {
                                                   >> 691           inside = true;
                                                   >> 692           latDisplasment = false;
                                                   >> 693           return ConvertTrueToGeom(tPathLength, currentMinimalStep);
                                                   >> 694         }
                                                   >> 695 
                                                   >> 696       if(firstStep || stepStatus == fGeomBoundary)
                                                   >> 697       {
618         rangeinit = currentRange;                 698         rangeinit = currentRange;
619         fr = facrange;                            699         fr = facrange;
                                                   >> 700         rangecut = geombig;
620         if(mass < masslimite)                     701         if(mass < masslimite)
621           {                                    << 702         {
622             if(lambda0 > lambdalimit) {        << 703           G4int index = 1;
623               fr *= (0.84+0.16*lambda0/lambdal << 704           if(charge > 0.) index = 2;
624             }                                  << 705           rangecut = couple->GetProductionCuts()->GetProductionCut(index);
                                                   >> 706           if(lambda0 > lambdalimit) {
                                                   >> 707             fr *= (0.84+0.16*lambda0/lambdalimit);
625           }                                       708           }
                                                   >> 709         }
                                                   >> 710 
626         //lower limit for tlimit                  711         //lower limit for tlimit
627         stepmin = ComputeStepmin();            << 712         G4double rat = currentKinEnergy/MeV;
628         tlimitmin = ComputeTlimitmin();        << 713         rat = 1.e-3/(rat*(10 + rat)) ;
                                                   >> 714         stepmin = lambda0*rat;
                                                   >> 715         tlimitmin = max(10*stepmin, tlimitminfix);
629       }                                           716       }
630       //step limit                                717       //step limit
631       tlimit = (currentRange > presafety) ?    << 718       tlimit = max(fr*rangeinit, facsafety*presafety);
632   std::max(fr*rangeinit, facsafety*presafety)  << 
633                                                   719 
634       //lower limit for tlimit                    720       //lower limit for tlimit
635       tlimit = std::max(tlimit, tlimitmin);    << 721       tlimit = max(tlimit, tlimitmin);
636                                                   722 
637       // condition for tPathLength from drr an    723       // condition for tPathLength from drr and finalr
638       if(currentRange > finalr) {                 724       if(currentRange > finalr) {
639         G4double tmax = drr*currentRange+         725         G4double tmax = drr*currentRange+
640                         finalr*(1.-drr)*(2.-fi    726                         finalr*(1.-drr)*(2.-finalr/currentRange);
641         tPathLength = std::min(tPathLength,tma << 727         tPathLength = min(tPathLength,tmax); 
642       }                                           728       }
643                                                   729 
644       // randomise if step determined by msc   << 730       // condition safety
645       tPathLength = (tlimit < tPathLength) ?   << 731       if(currentRange > rangecut) {
646         std::min(tPathLength, Randomizetlimit( << 732         if(firstStep) {
                                                   >> 733           tPathLength = min(tPathLength,facsafety*presafety);
                                                   >> 734         }
                                                   >> 735       else if(stepStatus != fGeomBoundary) {
                                                   >> 736         if(presafety > stepmin) {
                                                   >> 737           tPathLength = min(tPathLength,presafety);
                                                   >> 738           }
                                                   >> 739         }
                                                   >> 740       } 
                                                   >> 741 
                                                   >> 742       if(firstStep || stepStatus == fGeomBoundary)
                                                   >> 743       {
                                                   >> 744         G4double temptlimit = tlimit;
                                                   >> 745         if(temptlimit > tlimitmin)
                                                   >> 746         {
                                                   >> 747           do {
                                                   >> 748             temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.3*tlimit);
                                                   >> 749           } while ((temptlimit < tlimitmin) ||
                                                   >> 750                    (temptlimit > 2.*tlimit-tlimitmin));
                                                   >> 751         }
                                                   >> 752         else { temptlimit = tlimitmin; }
                                                   >> 753 
                                                   >> 754         tPathLength = min(tPathLength, temptlimit);
                                                   >> 755       }
                                                   >> 756       else { tPathLength = min(tPathLength, tlimit); }
647     }                                             757     }
648                                                   758 
649   // ----------------------------------------- << 759   // version similar to 7.1 (needed for some experiments)
650   // simple step limitation                    << 
651   else                                            760   else
652     {                                             761     {
653       if (stepStatus == fGeomBoundary)            762       if (stepStatus == fGeomBoundary)
                                                   >> 763   {
                                                   >> 764     if (currentRange > lambda0) { tlimit = facrange*currentRange; }
                                                   >> 765     else                        { tlimit = facrange*lambda0; }
                                                   >> 766 
                                                   >> 767           tlimit = max(tlimit, tlimitmin);
                                                   >> 768         }
                                                   >> 769 
                                                   >> 770       if(firstStep || stepStatus == fGeomBoundary)
                                                   >> 771       {
                                                   >> 772         G4double temptlimit = tlimit;
                                                   >> 773         if(temptlimit > tlimitmin)
654         {                                         774         {
655           tlimit = (currentRange > lambda0)    << 775           do {
656       ? facrange*currentRange : facrange*lambd << 776             temptlimit = G4RandGauss::shoot(rndmEngineMod,tlimit,0.3*tlimit); 
657           tlimit = std::max(tlimit, tlimitmin) << 777     } while ((temptlimit < tlimitmin) ||
                                                   >> 778        (temptlimit > 2.*tlimit-tlimitmin));
658         }                                         779         }
659       // randomise if step determined by msc   << 780         else { temptlimit = tlimitmin; }
660       tPathLength = (tlimit < tPathLength) ?   << 
661         std::min(tPathLength, Randomizetlimit( << 
662     }                                          << 
663                                                   781 
664   // ----------------------------------------- << 782         tPathLength = min(tPathLength, temptlimit);
665   firstStep = false;                           << 783       }
                                                   >> 784       else { tPathLength = min(tPathLength, tlimit); }
                                                   >> 785     }
                                                   >> 786   //G4cout << "tPathLength= " << tPathLength 
                                                   >> 787   //   << " currentMinimalStep= " << currentMinimalStep << G4endl;
666   return ConvertTrueToGeom(tPathLength, curren    788   return ConvertTrueToGeom(tPathLength, currentMinimalStep);
667 }                                                 789 }
668                                                   790 
669 //....oooOO0OOooo........oooOO0OOooo........oo    791 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
670                                                   792 
671 G4double G4UrbanMscModel::ComputeGeomPathLengt    793 G4double G4UrbanMscModel::ComputeGeomPathLength(G4double)
672 {                                                 794 {
                                                   >> 795   firstStep = false; 
673   lambdaeff = lambda0;                            796   lambdaeff = lambda0;
674   par1 = -1. ;                                    797   par1 = -1. ;  
675   par2 = par3 = 0. ;                              798   par2 = par3 = 0. ;  
676                                                   799 
677   // this correction needed to run MSC with eI    800   // this correction needed to run MSC with eIoni and eBrem inactivated
678   // and makes no harm for a normal run           801   // and makes no harm for a normal run
679   tPathLength = std::min(tPathLength,currentRa << 802   tPathLength = min(tPathLength,currentRange); 
680                                                   803 
681   //  do the true -> geom transformation          804   //  do the true -> geom transformation
682   zPathLength = tPathLength;                      805   zPathLength = tPathLength;
683                                                   806 
684   // z = t for very small tPathLength             807   // z = t for very small tPathLength
685   if(tPathLength < tlimitminfix2) return zPath    808   if(tPathLength < tlimitminfix2) return zPathLength;
686                                                   809 
                                                   >> 810   // VI: it is already checked
                                                   >> 811   // if(tPathLength > currentRange)
                                                   >> 812   //  tPathLength = currentRange ;
687   /*                                              813   /*
688   G4cout << "ComputeGeomPathLength: tpl= " <<     814   G4cout << "ComputeGeomPathLength: tpl= " <<  tPathLength
689          << " R= " << currentRange << " L0= "  << 815    << " R= " << currentRange << " L0= " << lambda0
690          << " E= " << currentKinEnergy << "  " << 816    << " E= " << currentKinEnergy << "  " 
691          << particle->GetParticleName() << G4e << 817    << particle->GetParticleName() << G4endl;
692   */                                              818   */
693   G4double tau = tPathLength/lambda0 ;            819   G4double tau = tPathLength/lambda0 ;
694                                                   820 
695   if ((tau <= tausmall) || insideskin) {          821   if ((tau <= tausmall) || insideskin) {
696     zPathLength = std::min(tPathLength, lambda << 822     zPathLength = min(tPathLength, lambda0); 
697                                                   823 
698   } else if (tPathLength < currentRange*dtrl)  << 824   } else  if (tPathLength < currentRange*dtrl) {
699     zPathLength = (tau < taulim) ? tPathLength << 825     if(tau < taulim) zPathLength = tPathLength*(1.-0.5*tau) ;
700       : lambda0*(1.-G4Exp(-tau));              << 826     else             zPathLength = lambda0*(1.-G4Exp(-tau));
701                                                   827 
702   } else if(currentKinEnergy < mass || tPathLe    828   } else if(currentKinEnergy < mass || tPathLength == currentRange)  {
703     par1 = 1./currentRange;                    << 829     par1 = 1./currentRange ;
704     par2 = currentRange/lambda0;               << 830     par2 = 1./(par1*lambda0) ;
705     par3 = 1.+par2;                            << 831     par3 = 1.+par2 ;
706     if(tPathLength < currentRange) {              832     if(tPathLength < currentRange) {
707       zPathLength =                               833       zPathLength = 
708         (1.-G4Exp(par3*G4Log(1.-tPathLength/cu << 834   (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3);
709     } else {                                      835     } else {
710       zPathLength = 1./(par1*par3);               836       zPathLength = 1./(par1*par3);
711     }                                             837     }
712                                                   838 
713   } else {                                        839   } else {
714     G4double rfin = std::max(currentRange-tPat << 840     G4double rfin = max(currentRange-tPathLength, 0.01*currentRange);
715     G4double T1 = GetEnergy(particle,rfin,coup    841     G4double T1 = GetEnergy(particle,rfin,couple);
716     G4double lambda1 = GetTransportMeanFreePat    842     G4double lambda1 = GetTransportMeanFreePath(particle,T1);
717                                                   843 
718     par1 = (lambda0-lambda1)/(lambda0*tPathLen    844     par1 = (lambda0-lambda1)/(lambda0*tPathLength);
719     //G4cout << "par1= " << par1 << " L1= " <<    845     //G4cout << "par1= " << par1 << " L1= " << lambda1 << G4endl;
720     par2 = 1./(par1*lambda0);                     846     par2 = 1./(par1*lambda0);
721     par3 = 1.+par2;                            << 847     par3 = 1.+par2 ;
722     zPathLength = (1.-G4Exp(par3*G4Log(lambda1    848     zPathLength = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3);
723   }                                               849   }
724                                                   850 
725   zPathLength = std::min(zPathLength, lambda0) << 851   zPathLength = min(zPathLength, lambda0);
726   //G4cout<< "zPathLength= "<< zPathLength<< "    852   //G4cout<< "zPathLength= "<< zPathLength<< " L0= " << lambda0 << G4endl;
727   return zPathLength;                             853   return zPathLength;
728 }                                                 854 }
729                                                   855 
730 //....oooOO0OOooo........oooOO0OOooo........oo    856 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
731                                                   857 
732 G4double G4UrbanMscModel::ComputeTrueStepLengt    858 G4double G4UrbanMscModel::ComputeTrueStepLength(G4double geomStepLength)
733 {                                                 859 {
734   // step defined other than transportation       860   // step defined other than transportation 
735   if(geomStepLength == zPathLength) {             861   if(geomStepLength == zPathLength) { 
736     //G4cout << "Urban::ComputeTrueLength: tPa    862     //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength 
737     //           << " step= " << geomStepLengt << 863     //     << " step= " << geomStepLength << " *** " << G4endl;
738     return tPathLength;                           864     return tPathLength; 
739   }                                               865   }
740                                                   866 
741   zPathLength = geomStepLength;                   867   zPathLength = geomStepLength;
742                                                   868 
743   // t = z for very small step                    869   // t = z for very small step
744   if(geomStepLength < tlimitminfix2) {            870   if(geomStepLength < tlimitminfix2) { 
745     tPathLength = geomStepLength;                 871     tPathLength = geomStepLength; 
746                                                   872   
747   // recalculation                                873   // recalculation
748   } else {                                        874   } else {
749                                                   875 
750     G4double tlength = geomStepLength;            876     G4double tlength = geomStepLength;
751     if((geomStepLength > lambda0*tausmall) &&     877     if((geomStepLength > lambda0*tausmall) && !insideskin) {
752                                                   878 
753       if(par1 <  0.) {                            879       if(par1 <  0.) {
754         tlength = -lambda0*G4Log(1.-geomStepLe << 880   tlength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
755       } else {                                    881       } else {
756         const G4double par4 = par1*par3;       << 882   if(par1*par3*geomStepLength < 1.) {
757         if(par4*geomStepLength < 1.) {         << 883     tlength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
758           tlength = (1.-G4Exp(G4Log(1.-par4*ge << 884   } else {
759         } else {                               << 885     tlength = currentRange;
760           tlength = currentRange;              << 886   }
761         }                                      << 
762       }                                           887       }
763                                                   888 
764       if(tlength < geomStepLength)   { tlength    889       if(tlength < geomStepLength)   { tlength = geomStepLength; }
765       else if(tlength > tPathLength) { tlength    890       else if(tlength > tPathLength) { tlength = tPathLength; }
766     }                                             891     }  
767     tPathLength = tlength;                        892     tPathLength = tlength; 
768   }                                               893   }
769   //G4cout << "Urban::ComputeTrueLength: tPath    894   //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength 
770   //         << " step= " << geomStepLength << << 895   //   << " step= " << geomStepLength << " &&& " << G4endl;
771                                                   896 
772   return tPathLength;                             897   return tPathLength;
773 }                                                 898 }
774                                                   899 
775 //....oooOO0OOooo........oooOO0OOooo........oo    900 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
776                                                   901 
777 G4ThreeVector&                                    902 G4ThreeVector& 
778 G4UrbanMscModel::SampleScattering(const G4Thre    903 G4UrbanMscModel::SampleScattering(const G4ThreeVector& oldDirection,
779                                   G4double /*s << 904           G4double /*safety*/)
780 {                                                 905 {
781   fDisplacement.set(0.0,0.0,0.0);                 906   fDisplacement.set(0.0,0.0,0.0);
782   if(tPathLength >= currentRange) { return fDi << 907   G4double kineticEnergy = currentKinEnergy;
783                                                << 
784   G4double kinEnergy = currentKinEnergy;       << 
785   if (tPathLength > currentRange*dtrl) {          908   if (tPathLength > currentRange*dtrl) {
786     kinEnergy = GetEnergy(particle,currentRang << 909     kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
787   } else if(tPathLength > currentRange*0.01) { << 910   } else {
788     kinEnergy -= tPathLength*GetDEDX(particle, << 911     kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
789                                      currentLo << 
790   }                                               912   }
791                                                   913 
792   if((tPathLength <= tlimitminfix) || (tPathLe << 914   if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
793      (kinEnergy <= CLHEP::eV)) { return fDispl << 915      (tPathLength < tausmall*lambda0)) { return fDisplacement; }
794                                                   916 
795   G4double cth = SampleCosineTheta(tPathLength << 917   G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
796                                                   918 
797   // protection against 'bad' cth values          919   // protection against 'bad' cth values
798   if(std::abs(cth) >= 1.0) { return fDisplacem << 920   if(std::fabs(cth) >= 1.0) { return fDisplacement; } 
799                                                   921 
800   G4double sth = std::sqrt((1.0 - cth)*(1.0 +  << 922   /*
801   G4double phi = CLHEP::twopi*rndmEngineMod->f << 923   if(cth < 1.0 - 1000*tPathLength/lambda0 && cth < 0.5 &&
802   G4ThreeVector newDirection(sth*std::cos(phi) << 924      kineticEnergy > 20*MeV) { 
803   newDirection.rotateUz(oldDirection);         << 925     G4cout << "### G4UrbanMscModel::SampleScattering for "
                                                   >> 926      << particle->GetParticleName()
                                                   >> 927      << " E(MeV)= " << kineticEnergy/MeV
                                                   >> 928      << " Step(mm)= " << tPathLength/mm
                                                   >> 929      << " in " << CurrentCouple()->GetMaterial()->GetName()
                                                   >> 930      << " CosTheta= " << cth << G4endl;
                                                   >> 931   }
                                                   >> 932   */
                                                   >> 933   G4double sth  = sqrt((1.0 - cth)*(1.0 + cth));
                                                   >> 934   G4double phi  = twopi*rndmEngineMod->flat(); 
                                                   >> 935   G4double dirx = sth*cos(phi);
                                                   >> 936   G4double diry = sth*sin(phi);
804                                                   937 
                                                   >> 938   G4ThreeVector newDirection(dirx,diry,cth);
                                                   >> 939   newDirection.rotateUz(oldDirection);
805   fParticleChange->ProposeMomentumDirection(ne    940   fParticleChange->ProposeMomentumDirection(newDirection);
806   /*                                              941   /*
807   G4cout << "G4UrbanMscModel::SampleSecondarie    942   G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
808          << " sinTheta= " << sth << " safety(m << 943    << " sinTheta= " << sth << " safety(mm)= " << safety
809          << " trueStep(mm)= " << tPathLength   << 944    << " trueStep(mm)= " << tPathLength
810          << " geomStep(mm)= " << zPathLength   << 945    << " geomStep(mm)= " << zPathLength
811          << G4endl;                            << 946    << G4endl;
812   */                                              947   */
813                                                   948 
                                                   >> 949   //if (latDisplasment && safety > tlimitminfix2 && currentTau >= tausmall &&
                                                   >> 950   //    !insideskin) {
814   if (latDisplasment && currentTau >= tausmall    951   if (latDisplasment && currentTau >= tausmall) {
815     if(dispAlg96) { SampleDisplacement(sth, ph << 952     //sample displacement r
816     else          { SampleDisplacementNew(cth, << 953     G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
817     fDisplacement.rotateUz(oldDirection);      << 954     G4double r = rmax*G4Exp(G4Log(rndmEngineMod->flat())*third);
                                                   >> 955 
                                                   >> 956     /*    
                                                   >> 957     G4cout << "G4UrbanMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
                                                   >> 958      << " sinTheta= " << sth << " r(mm)= " << r
                                                   >> 959            << " trueStep(mm)= " << tPathLength
                                                   >> 960            << " geomStep(mm)= " << zPathLength
                                                   >> 961            << G4endl;
                                                   >> 962     */
                                                   >> 963     if(r > 0.)
                                                   >> 964       {
                                                   >> 965         G4double latcorr = LatCorrelation();
                                                   >> 966         latcorr = min(latcorr, r);
                                                   >> 967 
                                                   >> 968         // sample direction of lateral displacement
                                                   >> 969         // compute it from the lateral correlation
                                                   >> 970         G4double Phi = 0.;
                                                   >> 971         if(std::abs(r*sth) < latcorr)
                                                   >> 972           Phi  = twopi*rndmEngineMod->flat();
                                                   >> 973         else
                                                   >> 974         {
                                                   >> 975     //G4cout << "latcorr= " << latcorr << "  r*sth= " << r*sth 
                                                   >> 976     //   << " ratio= " << latcorr/(r*sth) <<  G4endl;
                                                   >> 977           G4double psi = std::acos(latcorr/(r*sth));
                                                   >> 978           if(rndmEngineMod->flat() < 0.5)
                                                   >> 979             Phi = phi+psi;
                                                   >> 980           else
                                                   >> 981             Phi = phi-psi;
                                                   >> 982         }
                                                   >> 983 
                                                   >> 984         dirx = std::cos(Phi);
                                                   >> 985         diry = std::sin(Phi);
                                                   >> 986 
                                                   >> 987         fDisplacement.set(r*dirx,r*diry,0.0);
                                                   >> 988         fDisplacement.rotateUz(oldDirection);
                                                   >> 989       }
818   }                                               990   }
819   return fDisplacement;                           991   return fDisplacement;
820 }                                                 992 }
821                                                   993 
822 //....oooOO0OOooo........oooOO0OOooo........oo    994 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
823                                                   995 
824 G4double G4UrbanMscModel::SampleCosineTheta(G4    996 G4double G4UrbanMscModel::SampleCosineTheta(G4double trueStepLength,
825                                             G4 << 997               G4double KineticEnergy)
826 {                                                 998 {
827   G4double cth = 1.0;                          << 999   G4double cth = 1. ;
828   G4double tau = trueStepLength/lambda0;          1000   G4double tau = trueStepLength/lambda0;
                                                   >> 1001   currentTau   = tau;
                                                   >> 1002   lambdaeff    = lambda0;
829                                                   1003 
830   // mean tau value                            << 1004   Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
831   if(currentKinEnergy != kinEnergy) {          << 1005          couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
832     G4double lambda1 = GetTransportMeanFreePat << 1006 
833     if(std::abs(lambda1 - lambda0) > lambda0*0 << 1007   if(Zold != Zeff)  
834       tau = trueStepLength*G4Log(lambda0/lambd << 1008     UpdateCache();
835     }                                          << 1009 
                                                   >> 1010   G4double lambda1 = GetTransportMeanFreePath(particle,KineticEnergy);
                                                   >> 1011   if(std::fabs(lambda1 - lambda0) > lambda0*0.01 && lambda1 > 0.)
                                                   >> 1012   {
                                                   >> 1013     // mean tau value
                                                   >> 1014     tau = trueStepLength*G4Log(lambda0/lambda1)/(lambda0-lambda1);
836   }                                               1015   }
837                                                   1016 
838   currentTau = tau;                            << 1017   currentTau = tau ;
839   lambdaeff = trueStepLength/currentTau;          1018   lambdaeff = trueStepLength/currentTau;
840   currentRadLength = couple->GetMaterial()->Ge    1019   currentRadLength = couple->GetMaterial()->GetRadlen();
841                                                   1020 
842   if (tau >= taubig) { cth = -1.+2.*rndmEngine    1021   if (tau >= taubig) { cth = -1.+2.*rndmEngineMod->flat(); }
843   else if (tau >= tausmall) {                     1022   else if (tau >= tausmall) {
844     static const G4double numlim = 0.01;          1023     static const G4double numlim = 0.01;
845     static const G4double onethird = 1./3.;    << 1024     G4double xmeanth, x2meanth;
846     if(tau < numlim) {                            1025     if(tau < numlim) {
847       xmeanth = 1.0 - tau*(1.0 - 0.5*tau);        1026       xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
848       x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*one << 1027       x2meanth= 1.0 - tau*(5.0 - 6.25*tau)/3.;
849     } else {                                      1028     } else {
850       xmeanth = G4Exp(-tau);                      1029       xmeanth = G4Exp(-tau);
851       x2meanth = (1.+2.*G4Exp(-2.5*tau))*oneth << 1030       x2meanth = (1.+2.*G4Exp(-2.5*tau))/3.;
852     }                                             1031     }
853                                                   1032 
854     // too large step of low-energy particle      1033     // too large step of low-energy particle
855     G4double relloss = 1. - kinEnergy/currentK << 1034     G4double relloss = 1. - KineticEnergy/currentKinEnergy;
856     static const G4double rellossmax= 0.50;    << 
857     if(relloss > rellossmax) {                    1035     if(relloss > rellossmax) {
858       return SimpleScattering();               << 1036       return SimpleScattering(xmeanth,x2meanth);
859     }                                             1037     }
860     // is step extreme small ?                    1038     // is step extreme small ?
861     G4bool extremesmallstep = false;           << 1039     G4bool extremesmallstep = false ;
862     G4double tsmall = std::min(tlimitmin,lambd << 1040     G4double tsmall = tlimitmin;
863                                                << 1041     G4double theta0 = 0.;
864     G4double theta0;                           << 
865     if(trueStepLength > tsmall) {                 1042     if(trueStepLength > tsmall) {
866       theta0 = ComputeTheta0(trueStepLength,ki << 1043        theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
867     } else {                                      1044     } else {
868       theta0 = std::sqrt(trueStepLength/tsmall << 1045       G4double rate = trueStepLength/tsmall ;
869   *ComputeTheta0(tsmall,kinEnergy);            << 1046       if(rndmEngineMod->flat() < rate) {
870       extremesmallstep = true;                 << 1047   theta0 = ComputeTheta0(tsmall,KineticEnergy);
                                                   >> 1048   extremesmallstep = true ;
                                                   >> 1049       }
871     }                                             1050     }
872                                                << 1051     //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max 
873     static const G4double onesixth = 1./6.;    << 1052     //       << "  sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
874     static const G4double one12th = 1./12.;    << 
875     static const G4double theta0max = CLHEP::p << 
876                                                   1053 
877     // protection for very small angles           1054     // protection for very small angles
878     G4double theta2 = theta0*theta0;              1055     G4double theta2 = theta0*theta0;
879                                                   1056 
880     if(theta2 < tausmall) { return cth; }         1057     if(theta2 < tausmall) { return cth; }
881     if(theta0 > theta0max) { return SimpleScat << 1058     
                                                   >> 1059     if(theta0 > theta0max) {
                                                   >> 1060       return SimpleScattering(xmeanth,x2meanth);
                                                   >> 1061     }
882                                                   1062 
883     G4double x = theta2*(1.0 - theta2*one12th) << 1063     G4double x = theta2*(1.0 - theta2/12.);
884     if(theta2 > numlim) {                         1064     if(theta2 > numlim) {
885       G4double sth = 2*std::sin(0.5*theta0);   << 1065       G4double sth = 2*sin(0.5*theta0);
886       x = sth*sth;                                1066       x = sth*sth;
887     }                                             1067     }
888                                                   1068 
889     // parameter for tail                         1069     // parameter for tail
890     G4double ltau = G4Log(tau);                << 1070     G4double ltau= G4Log(tau);
891     G4double u = !extremesmallstep ? G4Exp(lta << 1071     G4double u   = G4Exp(ltau/6.);
892       : G4Exp(G4Log(tsmall/lambda0)*onesixth); << 1072     if(extremesmallstep)  u = G4Exp(G4Log(tsmall/lambda0)/6.);
893                                                << 
894     G4double xx  = G4Log(lambdaeff/currentRadL    1073     G4double xx  = G4Log(lambdaeff/currentRadLength);
895     G4double xsi = msc[idx]->coeffc1 +         << 1074     G4double xsi = coeffc1+u*(coeffc2+coeffc3*u)+coeffc4*xx;
896       u*(msc[idx]->coeffc2+msc[idx]->coeffc3*u << 
897                                                   1075 
898     // tail should not be too big                 1076     // tail should not be too big
899     xsi = std::max(xsi, 1.9);                  << 1077     if(xsi < 1.9) { 
900       /*                                          1078       /*
901       if(KineticEnergy > 20*MeV && xsi < 1.6)     1079       if(KineticEnergy > 20*MeV && xsi < 1.6) {
902         G4cout << "G4UrbanMscModel::SampleCosi << 1080   G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= " 
903                << KineticEnergy/GeV            << 1081          << KineticEnergy/GeV 
904                << " !!** c= " << xsi           << 1082          << " !!** c= " << xsi
905                << " **!! length(mm)= " << true << 1083          << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff 
906                << " " << couple->GetMaterial() << 1084          << " " << couple->GetMaterial()->GetName()
907                << " tau= " << tau << G4endl;   << 1085          << " tau= " << tau << G4endl;
908       }                                           1086       }
909       */                                          1087       */
                                                   >> 1088       xsi = 1.9; 
                                                   >> 1089     }
910                                                   1090 
911     G4double c = xsi;                             1091     G4double c = xsi;
912                                                   1092 
913     if(std::abs(c-3.) < 0.001)      { c = 3.00 << 1093     if(fabs(c-3.) < 0.001)      { c = 3.001; }
914     else if(std::abs(c-2.) < 0.001) { c = 2.00 << 1094     else if(fabs(c-2.) < 0.001) { c = 2.001; }
915                                                   1095 
916     G4double c1 = c-1.;                           1096     G4double c1 = c-1.;
                                                   >> 1097 
917     G4double ea = G4Exp(-xsi);                    1098     G4double ea = G4Exp(-xsi);
918     G4double eaa = 1.-ea ;                        1099     G4double eaa = 1.-ea ;
919     G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/ea    1100     G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
920     G4double x0 = 1. - xsi*x;                     1101     G4double x0 = 1. - xsi*x;
921                                                   1102 
922     // G4cout << " xmean1= " << xmean1 << "  x    1103     // G4cout << " xmean1= " << xmean1 << "  xmeanth= " << xmeanth << G4endl;
923                                                   1104 
924     if(xmean1 <= 0.999*xmeanth) { return Simpl << 1105     if(xmean1 <= 0.999*xmeanth) {
925                                                << 1106       return SimpleScattering(xmeanth,x2meanth);
                                                   >> 1107     }
926     //from continuity of derivatives              1108     //from continuity of derivatives
927     G4double b = 1.+(c-xsi)*x;                    1109     G4double b = 1.+(c-xsi)*x;
928                                                   1110 
929     G4double b1 = b+1.;                           1111     G4double b1 = b+1.;
930     G4double bx = c*x;                            1112     G4double bx = c*x;
931                                                   1113 
932     G4double eb1 = G4Exp(G4Log(b1)*c1);           1114     G4double eb1 = G4Exp(G4Log(b1)*c1);
933     G4double ebx = G4Exp(G4Log(bx)*c1);           1115     G4double ebx = G4Exp(G4Log(bx)*c1);
934     G4double d = ebx/eb1;                         1116     G4double d = ebx/eb1;
935                                                   1117 
936     G4double xmean2 = (x0 + d - (bx - b1*d)/(c    1118     G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
937                                                   1119 
938     G4double f1x0 = ea/eaa;                       1120     G4double f1x0 = ea/eaa;
939     G4double f2x0 = c1/(c*(1. - d));              1121     G4double f2x0 = c1/(c*(1. - d));
940     G4double prob = f2x0/(f1x0+f2x0);             1122     G4double prob = f2x0/(f1x0+f2x0);
941                                                   1123 
942     G4double qprob = xmeanth/(prob*xmean1+(1.-    1124     G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
943                                                   1125 
944     // sampling of costheta                       1126     // sampling of costheta
945     //G4cout << "c= " << c << " qprob= " << qp    1127     //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
946     // << " c1= " << c1 << " b1= " << b1 << "     1128     // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
947     //             << G4endl;                  << 1129     //       << G4endl;
948     rndmEngineMod->flatArray(2, rndmarray);    << 1130     if(rndmEngineMod->flat() < qprob)
949     if(rndmarray[0] < qprob)                   << 
950     {                                             1131     {
951       G4double var = 0;                           1132       G4double var = 0;
952       if(rndmarray[1] < prob) {                << 1133       if(rndmEngineMod->flat() < prob) {
953         cth = 1.+G4Log(ea+rndmEngineMod->flat(    1134         cth = 1.+G4Log(ea+rndmEngineMod->flat()*eaa)*x;
954       } else {                                    1135       } else {
955         var = (1.0 - d)*rndmEngineMod->flat();    1136         var = (1.0 - d)*rndmEngineMod->flat();
956         if(var < numlim*d) {                      1137         if(var < numlim*d) {
957           var /= (d*c1);                          1138           var /= (d*c1); 
958           cth = -1.0 + var*(1.0 - 0.5*var*c)*(    1139           cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
959         } else {                                  1140         } else {
960           cth = 1. + x*(c - xsi - c*G4Exp(-G4L    1141           cth = 1. + x*(c - xsi - c*G4Exp(-G4Log(var + d)/c1));
961         }                                         1142         }
962       }                                           1143       } 
963     } else {                                   << 1144       /*
964       cth = -1.+2.*rndmarray[1];               << 1145       if(KineticEnergy > 5*GeV && cth < 0.9) {
965     }                                          << 1146         G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= " 
966   }                                            << 1147                << KineticEnergy/GeV 
967   return cth;                                  << 1148          << " 1-cosT= " << 1 - cth
968 }                                              << 1149                << " length(mm)= " << trueStepLength << " Zeff= " << Zeff 
969                                                << 1150                << " tau= " << tau
970 //....oooOO0OOooo........oooOO0OOooo........oo << 1151          << " prob= " << prob << " var= " << var << G4endl;
971                                                << 1152   G4cout << "  c= " << c << " qprob= " << qprob << " eb1= " << eb1
972 G4double G4UrbanMscModel::ComputeTheta0(G4doub << 1153          << " ebx= " << ebx
973                                         G4doub << 1154                << " c1= " << c1 << " b= " << b << " b1= " << b1 
974 {                                              << 1155                << " bx= " << bx << " d= " << d
975   // for all particles take the width of the c << 1156                << " ea= " << ea << " eaa= " << eaa << G4endl;
976   //  from a  parametrization similar to the H << 1157       }
977   // ( Highland formula: Particle Physics Book << 1158       */
978   G4double invbetacp = (kinEnergy+mass)/(kinEn << 
979   if(currentKinEnergy != kinEnergy) {          << 
980     invbetacp = std::sqrt(invbetacp*(currentKi << 
981         (currentKinEnergy*(currentKinEnergy+2. << 
982   }                                            << 
983   G4double y = trueStepLength/currentRadLength << 
984                                                << 
985   if(fPosiCorrection && particle == positron)  << 
986   {                                            << 
987     static const G4double xl= 0.6;             << 
988     static const G4double xh= 0.9;             << 
989     static const G4double e = 113.0;           << 
990     G4double corr;                             << 
991                                                << 
992     G4double tau = std::sqrt(currentKinEnergy* << 
993     G4double x = std::sqrt(tau*(tau+2.)/((tau+ << 
994     G4double a = msc[idx]->posa;               << 
995     G4double b = msc[idx]->posb;               << 
996     G4double c = msc[idx]->posc;               << 
997     G4double d = msc[idx]->posd;               << 
998     if(x < xl) {                               << 
999       corr = a*(1.-G4Exp(-b*x));               << 
1000     } else if(x > xh) {                       << 
1001       corr = c+d*G4Exp(e*(x-1.));             << 
1002     } else {                                  << 
1003       G4double yl = a*(1.-G4Exp(-b*xl));      << 
1004       G4double yh = c+d*G4Exp(e*(xh-1.));     << 
1005       G4double y0 = (yh-yl)/(xh-xl);          << 
1006       G4double y1 = yl-y0*xl;                 << 
1007       corr = y0*x+y1;                         << 
1008     }                                         << 
1009     //======================================= << 
1010     y *= corr*msc[idx]->pose;                 << 
1011   }                                           << 
1012                                               << 
1013   static const G4double c_highland = 13.6*CLH << 
1014   G4double theta0 = c_highland*std::abs(charg << 
1015                                               << 
1016   // correction factor from e- scattering dat << 
1017   theta0 *= (msc[idx]->coeffth1+msc[idx]->coe << 
1018   return theta0;                              << 
1019 }                                             << 
1020                                               << 
1021 //....oooOO0OOooo........oooOO0OOooo........o << 
1022                                               << 
1023 void G4UrbanMscModel::SampleDisplacement(G4do << 
1024 {                                             << 
1025   // simple and fast sampling                 << 
1026   // based on single scattering results       << 
1027   // u = r/rmax : mean value                  << 
1028                                               << 
1029   G4double rmax = std::sqrt((tPathLength-zPat << 
1030   if(rmax > 0.)                               << 
1031   {                                           << 
1032     G4double r = 0.73*rmax;                   << 
1033                                               << 
1034     // simple distribution for v=Phi-phi=psi  << 
1035     // beta determined from the requirement t << 
1036     // the same mean value than that obtained << 
1037                                               << 
1038     static const G4double cbeta  = 2.160;     << 
1039     static const G4double cbeta1 = 1. - G4Exp << 
1040     rndmEngineMod->flatArray(2, rndmarray);   << 
1041     G4double psi = -G4Log(1. - rndmarray[0]*c << 
1042     G4double Phi = (rndmarray[1] < 0.5) ? phi << 
1043     fDisplacement.set(r*std::cos(Phi),r*std:: << 
1044   }                                           << 
1045 }                                             << 
1046                                               << 
1047 //....oooOO0OOooo........oooOO0OOooo........o << 
1048                                               << 
1049 void G4UrbanMscModel::SampleDisplacementNew(G << 
1050 {                                             << 
1051   // simple and fast sampling                 << 
1052   // based on single scattering results       << 
1053   // u = (r/rmax)**2 : distribution from ss s << 
1054   const G4double eps = 1.e-3;                 << 
1055   const G4double rmax =                       << 
1056     std::sqrt((tPathLength-zPathLength)*(tPat << 
1057                                               << 
1058   if(rmax > 0.)                               << 
1059   {                                           << 
1060     const G4double x0 = 0.73 ;                << 
1061     const G4double alpha = G4Log(7.33)/x0 ;   << 
1062     const G4double a1 = 1.-x0 ;               << 
1063     const G4double a2 = 1.-G4Exp(-alpha*x0) ; << 
1064     const G4double a3 = G4Exp(alpha*x0)-1. ;  << 
1065     const G4double w1 = 2.*a2/(alpha*a1+2.*a2 << 
1066                                               << 
1067     G4double r, sqx;                          << 
1068     if (rmax/currentRange < eps)              << 
1069     {                                         << 
1070       r = 0.73*rmax ;                         << 
1071       sqx = 1.;                               << 
1072     }                                            1159     }
1073     else                                      << 1160     else {
1074     {                                         << 1161       cth = -1.+2.*rndmEngineMod->flat();
1075       rndmEngineMod->flatArray(2,rndmarray);  << 1162       /*
1076       const G4double x = (rndmarray[0] < w1)  << 1163       if(KineticEnergy > 5*GeV) {
1077         1. - a1*std::sqrt(1.-rndmarray[1]);   << 1164         G4cout << "G4UrbanMscModel::SampleCosineTheta: E(GeV)= " 
1078                                               << 1165                << KineticEnergy/GeV 
1079       sqx = std::sqrt(x);                     << 1166                << " length(mm)= " << trueStepLength << " Zeff= " << Zeff 
1080       r = sqx*rmax;                           << 1167                << " qprob= " << qprob << G4endl;
                                                   >> 1168       }
                                                   >> 1169       */
1081     }                                            1170     }
1082     // Gaussian distribution for Phi-phi=psi  << 
1083     const G4double sigma = 0.1+0.9*sqx;       << 
1084     const G4double psi = G4RandGauss::shoot(0 << 
1085     const G4double Phi = phi+psi;             << 
1086     fDisplacement.set(r*std::cos(Phi), r*std: << 
1087   }                                              1171   }
                                                   >> 1172   return cth ;
1088 }                                                1173 }
1089                                                  1174 
1090 //....oooOO0OOooo........oooOO0OOooo........o    1175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1091                                                  1176 
1092 void G4UrbanMscModel::InitialiseModelCache()  << 
1093 {                                             << 
1094   // it is assumed, that for the second run o << 
1095   // of a new G4MaterialCutsCouple is possibl << 
1096   auto theCoupleTable = G4ProductionCutsTable << 
1097   std::size_t numOfCouples = theCoupleTable-> << 
1098   if(numOfCouples != msc.size()) { msc.resize << 
1099                                               << 
1100   for(G4int j=0; j<(G4int)numOfCouples; ++j)  << 
1101     auto aCouple = theCoupleTable->GetMateria << 
1102                                                  1177 
1103     // new couple                             << 
1104     msc[j] = new mscData();                   << 
1105     G4double Zeff = aCouple->GetMaterial()->G << 
1106     G4double sqrz = std::sqrt(Zeff);          << 
1107     msc[j]->sqrtZ = sqrz;                     << 
1108     // parameterisation of step limitation    << 
1109     msc[j]->factmin = dispAlg96 ? 0.001 : 0.0 << 
1110     G4double lnZ = G4Log(Zeff);               << 
1111     // correction in theta0 formula           << 
1112     G4double w = G4Exp(lnZ/6.);               << 
1113     G4double facz = 0.990395+w*(-0.168386+w*0 << 
1114     msc[j]->coeffth1 = facz*(1. - 8.7780e-2/Z << 
1115     msc[j]->coeffth2 = facz*(4.0780e-2 + 1.73 << 
1116                                               << 
1117     // tail parameters                        << 
1118     G4double Z13 = w*w;                       << 
1119     msc[j]->coeffc1 = 2.3785    - Z13*(4.1981 << 
1120     msc[j]->coeffc2 = 4.7526e-1 + Z13*(1.7694 << 
1121     msc[j]->coeffc3 = 2.3683e-1 - Z13*(1.8111 << 
1122     msc[j]->coeffc4 = 1.7888e-2 + Z13*(1.9659 << 
1123                                               << 
1124     msc[j]->Z23 = Z13*Z13;                    << 
1125                                               << 
1126     msc[j]->stepmina = 27.725/(1.+0.203*Zeff) << 
1127     msc[j]->stepminb =  6.152/(1.+0.111*Zeff) << 
1128                                               << 
1129     // 21.07.2020                             << 
1130     msc[j]->doverra = 9.6280e-1 - 8.4848e-2*m << 
1131                                               << 
1132     // 06.10.2020                             << 
1133     // msc[j]->doverra = 7.7024e-1 - 6.7878e- << 
1134     msc[j]->doverrb = 1.15 - 9.76e-4*Zeff;    << 
1135                                               << 
1136     // corrections for e+                     << 
1137     msc[j]->posa = 0.994-4.08e-3*Zeff;        << 
1138     msc[j]->posb = 7.16+(52.6+365./Zeff)/Zeff << 
1139     msc[j]->posc = 1.000-4.47e-3*Zeff;        << 
1140     msc[j]->posd = 1.21e-3*Zeff;              << 
1141     msc[j]->pose = 1.+Zeff*(1.84035e-4*Zeff-1 << 
1142   }                                           << 
1143 }                                             << 
1144                                                  1178 
1145 //....oooOO0OOooo........oooOO0OOooo........o << 
1146                                                  1179