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