| 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: $
>> 27 // GEANT4 tag $Name: $
>> 28 //
26 // ------------------------------------------- 29 // -------------------------------------------------------------------
>> 30 //
>> 31 // GEANT4 Class file
>> 32 //
27 // 33 //
28 // File name: G4UrbanAdjointMscModel 34 // File name: G4UrbanAdjointMscModel
29 // 35 //
30 // Author: Laszlo Urban 36 // Author: Laszlo Urban
31 // 37 //
32 // Creation date: 19.02.2013 38 // Creation date: 19.02.2013
33 // 39 //
34 // Created from G4UrbanAdjointMscModel96 40 // Created from G4UrbanAdjointMscModel96
35 // 41 //
>> 42 // New parametrization for theta0
>> 43 // Correction for very small step length
>> 44 //
>> 45 // Class Description:
>> 46 //
36 // Implementation of the model of multiple sca 47 // Implementation of the model of multiple scattering based on
37 // H.W.Lewis Phys Rev 78 (1950) 526 and others 48 // H.W.Lewis Phys Rev 78 (1950) 526 and others
38 // In its present form the model can be used << 49
39 // of the e-/e+ multiple scattering <<
40 // ------------------------------------------- 50 // -------------------------------------------------------------------
>> 51 // In its present form the model can be used for simulation
>> 52 // of the e-/e+ multiple scattering
>> 53 //
41 54
42 #include "G4UrbanAdjointMscModel.hh" <<
43 55
44 #include "globals.hh" << 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 58
>> 59 #include "G4UrbanAdjointMscModel.hh"
>> 60 #include "G4PhysicalConstants.hh"
>> 61 #include "G4SystemOfUnits.hh"
>> 62 #include "Randomize.hh"
45 #include "G4Electron.hh" 63 #include "G4Electron.hh"
46 #include "G4Exp.hh" << 64 #include "G4Positron.hh"
47 #include "G4Log.hh" <<
48 #include "G4LossTableManager.hh" 65 #include "G4LossTableManager.hh"
49 #include "G4ParticleChangeForMSC.hh" 66 #include "G4ParticleChangeForMSC.hh"
50 #include "G4PhysicalConstants.hh" << 67
51 #include "G4Poisson.hh" 68 #include "G4Poisson.hh"
52 #include "G4Positron.hh" <<
53 #include "G4Pow.hh" 69 #include "G4Pow.hh"
54 #include "G4SystemOfUnits.hh" << 70 #include "globals.hh"
55 #include "Randomize.hh" << 71 #include "G4Log.hh"
>> 72 #include "G4Exp.hh"
56 73
57 //....oooOO0OOooo........oooOO0OOooo........oo 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
58 75
59 using namespace std; 76 using namespace std;
60 77
61 //....oooOO0OOooo........oooOO0OOooo........oo 78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
62 79
63 G4UrbanAdjointMscModel::G4UrbanAdjointMscModel 80 G4UrbanAdjointMscModel::G4UrbanAdjointMscModel(const G4String& nam)
64 : G4VMscModel(nam) 81 : G4VMscModel(nam)
65 { 82 {
66 masslimite = 0.6 * MeV; << 83 masslimite = 0.6*MeV;
67 lambdalimit = 1. * mm; << 84 lambdalimit = 1.*mm;
68 fr = 0.02; 85 fr = 0.02;
69 taubig = 8.0; 86 taubig = 8.0;
70 tausmall = 1.e-16; 87 tausmall = 1.e-16;
71 taulim = 1.e-6; 88 taulim = 1.e-6;
72 currentTau = taulim; 89 currentTau = taulim;
73 tlimitminfix = 0.01 * nm; << 90 tlimitminfix = 0.01*nm;
74 tlimitminfix2 = 1. * nm; << 91 tlimitminfix2 = 1.*nm;
75 stepmin = tlimitminfix; 92 stepmin = tlimitminfix;
76 smallstep = 1.e10; 93 smallstep = 1.e10;
77 currentRange = 0.; << 94 currentRange = 0. ;
78 rangeinit = 0.; 95 rangeinit = 0.;
79 tlimit = 1.e10 * mm; << 96 tlimit = 1.e10*mm;
80 tlimitmin = 10. * tlimitminfix; << 97 tlimitmin = 10.*tlimitminfix;
81 tgeom = 1.e50 * mm; << 98 tgeom = 1.e50*mm;
82 geombig = 1.e50 * mm; << 99 geombig = 1.e50*mm;
83 geommin = 1.e-3 * mm; << 100 geommin = 1.e-3*mm;
84 geomlimit = geombig; 101 geomlimit = geombig;
85 presafety = 0. * mm; << 102 presafety = 0.*mm;
86 103
87 facsafety = 0.6; << 104 facsafety = 0.6;
88 105
89 Zold = 0.; << 106 Zold = 0.;
90 Zeff = 1.; << 107 Zeff = 1.;
91 Z2 = 1.; << 108 Z2 = 1.;
92 Z23 = 1.; << 109 Z23 = 1.;
93 lnZ = 0.; << 110 lnZ = 0.;
94 coeffth1 = 0.; << 111 coeffth1 = 0.;
95 coeffth2 = 0.; << 112 coeffth2 = 0.;
96 coeffc1 = 0.; << 113 coeffc1 = 0.;
97 coeffc2 = 0.; << 114 coeffc2 = 0.;
98 coeffc3 = 0.; << 115 coeffc3 = 0.;
99 coeffc4 = 0.; << 116 coeffc4 = 0.;
100 particle = nullptr; << 117 particle = 0;
101 118
102 positron = G4Positron::Positron(); 119 positron = G4Positron::Positron();
103 theManager = G4LossTableManager::Instance << 120 theManager = G4LossTableManager::Instance();
104 rndmEngineMod = G4Random::getTheEngine(); 121 rndmEngineMod = G4Random::getTheEngine();
105 122
106 firstStep = true; << 123 firstStep = true;
107 insideskin = false; << 124 insideskin = false;
108 latDisplasmentbackup = false; 125 latDisplasmentbackup = false;
109 displacementFlag = true; << 126 displacementFlag = true;
110 127
111 rangecut = geombig; 128 rangecut = geombig;
112 drr = 0.35; << 129 drr = 0.35 ;
113 finalr = 10. * um; << 130 finalr = 10.*um ;
114 131
115 skindepth = skin * stepmin; << 132 skindepth = skin*stepmin;
116 133
117 mass = proton_mass_c2; << 134 mass = proton_mass_c2;
118 charge = ChargeSquare = 1.0; 135 charge = ChargeSquare = 1.0;
119 currentKinEnergy = currentRadLength = lambda << 136 currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
120 zPathLength = par1 = par2 = par3 = 0.; << 137 = zPathLength = par1 = par2 = par3 = 0;
121 138
122 currentMaterialIndex = -1; 139 currentMaterialIndex = -1;
123 fParticleChange = nullptr; << 140 fParticleChange = nullptr;
124 couple = nullptr; << 141 couple = nullptr;
125 } 142 }
126 143
127 //....oooOO0OOooo........oooOO0OOooo........oo 144 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
128 G4UrbanAdjointMscModel::~G4UrbanAdjointMscMode << 145
>> 146 G4UrbanAdjointMscModel::~G4UrbanAdjointMscModel()
>> 147 {}
129 148
130 //....oooOO0OOooo........oooOO0OOooo........oo 149 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 150
131 void G4UrbanAdjointMscModel::Initialise(const 151 void G4UrbanAdjointMscModel::Initialise(const G4ParticleDefinition* p,
132 const << 152 const G4DataVector&)
133 { 153 {
134 const G4ParticleDefinition* p1 = p; << 154 const G4ParticleDefinition* p1 =p;
135 155
136 if(p->GetParticleName() == "adj_e-") << 156 if (p->GetParticleName() =="adj_e-") p1= G4Electron::Electron();
137 p1 = G4Electron::Electron(); <<
138 // set values of some data members 157 // set values of some data members
139 SetParticle(p1); 158 SetParticle(p1);
140 << 159 /*
>> 160 if(p->GetPDGMass() > MeV) {
>> 161 G4cout << "### WARNING: G4UrbanAdjointMscModel model is used for "
>> 162 << p->GetParticleName() << " !!! " << G4endl;
>> 163 G4cout << "### This model should be used only for e+-"
>> 164 << G4endl;
>> 165 }
>> 166 */
141 fParticleChange = GetParticleChangeForMSC(p1 167 fParticleChange = GetParticleChangeForMSC(p1);
142 168
143 latDisplasmentbackup = latDisplasment; 169 latDisplasmentbackup = latDisplasment;
>> 170
>> 171 //G4cout << "### G4UrbanAdjointMscModel::Initialise done!" << G4endl;
144 } 172 }
145 173
146 //....oooOO0OOooo........oooOO0OOooo........oo 174 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
147 175
148 G4double G4UrbanAdjointMscModel::ComputeCrossS 176 G4double G4UrbanAdjointMscModel::ComputeCrossSectionPerAtom(
149 const G4ParticleDefinition* part, G4double K << 177 const G4ParticleDefinition* part,
150 G4double AtomicNumber, G4double, G4double, G << 178 G4double KineticEnergy,
>> 179 G4double AtomicNumber,G4double,
>> 180 G4double, G4double)
151 { 181 {
152 static constexpr G4double epsmin = 1.e-4; << 182 static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
153 static constexpr G4double epsmax = 1.e10; <<
154 183
155 static constexpr G4double Zdat[15] = { 4., << 184 static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38.,47.,
156 47., << 185 50., 56., 64., 74., 79., 82. };
157 186
158 // corr. factors for e-/e+ lambda for T <= T 187 // corr. factors for e-/e+ lambda for T <= Tlim
159 static constexpr G4double celectron[15][22] << 188 static const G4double celectron[15][22] =
160 { 1.125, 1.072, 1.051, 1.047, 1.047, 1.050 << 189 {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
161 1.054, 1.057, 1.062, 1.069, 1.075, 1.090 << 190 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
162 1.112, 1.108, 1.100, 1.093, 1.089, 1.087 << 191 1.112,1.108,1.100,1.093,1.089,1.087 },
163 { 1.408, 1.246, 1.143, 1.096, 1.077, 1.059 << 192 {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
164 1.052, 1.053, 1.058, 1.065, 1.072, 1.087 << 193 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
165 1.109, 1.105, 1.097, 1.090, 1.086, 1.082 << 194 1.109,1.105,1.097,1.090,1.086,1.082 },
166 { 2.833, 2.268, 1.861, 1.612, 1.486, 1.309 << 195 {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
167 1.136, 1.114, 1.106, 1.106, 1.109, 1.119 << 196 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
168 1.131, 1.124, 1.113, 1.104, 1.099, 1.098 << 197 1.131,1.124,1.113,1.104,1.099,1.098 },
169 { 3.879, 3.016, 2.380, 2.007, 1.818, 1.535 << 198 {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
170 1.190, 1.133, 1.107, 1.099, 1.098, 1.103 << 199 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
171 1.112, 1.105, 1.096, 1.089, 1.085, 1.098 << 200 1.112,1.105,1.096,1.089,1.085,1.098 },
172 { 6.937, 4.330, 2.886, 2.256, 1.987, 1.628 << 201 {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
173 1.203, 1.122, 1.080, 1.065, 1.061, 1.063 << 202 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
174 1.073, 1.070, 1.064, 1.059, 1.056, 1.056 << 203 1.073,1.070,1.064,1.059,1.056,1.056 },
175 { 9.616, 5.708, 3.424, 2.551, 2.204, 1.762 << 204 {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
176 1.256, 1.155, 1.099, 1.077, 1.070, 1.068 << 205 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
177 1.074, 1.070, 1.063, 1.059, 1.056, 1.052 << 206 1.074,1.070,1.063,1.059,1.056,1.052 },
178 { 11.72, 6.364, 3.811, 2.806, 2.401, 1.884 << 207 {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
179 1.300, 1.180, 1.112, 1.082, 1.073, 1.066 << 208 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
180 1.068, 1.064, 1.059, 1.054, 1.051, 1.050 << 209 1.068,1.064,1.059,1.054,1.051,1.050 },
181 { 18.08, 8.601, 4.569, 3.183, 2.662, 2.025 << 210 {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
182 1.339, 1.195, 1.108, 1.068, 1.053, 1.040 << 211 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
183 1.039, 1.037, 1.034, 1.031, 1.030, 1.036 << 212 1.039,1.037,1.034,1.031,1.030,1.036 },
184 { 18.22, 10.48, 5.333, 3.713, 3.115, 2.367 << 213 {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
185 1.498, 1.301, 1.171, 1.105, 1.077, 1.048 << 214 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
186 1.031, 1.028, 1.024, 1.022, 1.021, 1.024 << 215 1.031,1.028,1.024,1.022,1.021,1.024 },
187 { 14.14, 10.65, 5.710, 3.929, 3.266, 2.453 << 216 {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
188 1.528, 1.319, 1.178, 1.106, 1.075, 1.040 << 217 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
189 1.020, 1.017, 1.015, 1.013, 1.013, 1.020 << 218 1.020,1.017,1.015,1.013,1.013,1.020 },
190 { 14.11, 11.73, 6.312, 4.240, 3.478, 2.566 << 219 {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
191 1.569, 1.342, 1.186, 1.102, 1.065, 1.022 << 220 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
192 0.995, 0.993, 0.993, 0.993, 0.993, 1.011 << 221 0.995,0.993,0.993,0.993,0.993,1.011 },
193 { 22.76, 20.01, 8.835, 5.287, 4.144, 2.901 << 222 {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
194 1.677, 1.410, 1.224, 1.121, 1.073, 1.014 << 223 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
195 0.974, 0.972, 0.973, 0.974, 0.975, 0.987 << 224 0.974,0.972,0.973,0.974,0.975,0.987 },
196 { 50.77, 40.85, 14.13, 7.184, 5.284, 3.435 << 225 {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
197 1.837, 1.512, 1.283, 1.153, 1.091, 1.010 << 226 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
198 0.950, 0.947, 0.949, 0.952, 0.954, 0.963 << 227 0.950,0.947,0.949,0.952,0.954,0.963 },
199 { 65.87, 59.06, 15.87, 7.570, 5.567, 3.650 << 228 {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
200 1.939, 1.579, 1.325, 1.178, 1.108, 1.014 << 229 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
201 0.941, 0.938, 0.940, 0.944, 0.946, 0.954 << 230 0.941,0.938,0.940,0.944,0.946,0.954 },
202 { 55.60, 47.34, 15.92, 7.810, 5.755, 3.767 << 231 {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
203 1.985, 1.609, 1.343, 1.188, 1.113, 1.013 << 232 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
204 0.933, 0.930, 0.933, 0.936, 0.939, 0.949 << 233 0.933,0.930,0.933,0.936,0.939,0.949 }};
205 }; << 234
206 << 235 static const G4double cpositron[15][22] = {
207 static constexpr G4double cpositron[15][22] << 236 {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
208 { 2.589, 2.044, 1.658, 1.446, 1.347, 1.217 << 237 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
209 1.097, 1.083, 1.080, 1.086, 1.092, 1.108 << 238 1.131,1.126,1.117,1.108,1.103,1.100 },
210 1.131, 1.126, 1.117, 1.108, 1.103, 1.100 << 239 {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
211 { 3.904, 2.794, 2.079, 1.710, 1.543, 1.325 << 240 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
212 1.122, 1.096, 1.089, 1.092, 1.098, 1.114 << 241 1.138,1.132,1.122,1.113,1.108,1.102 },
213 1.138, 1.132, 1.122, 1.113, 1.108, 1.102 << 242 {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
214 { 7.970, 6.080, 4.442, 3.398, 2.872, 2.127 << 243 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
215 1.357, 1.246, 1.194, 1.179, 1.178, 1.188 << 244 1.203,1.190,1.173,1.159,1.151,1.145 },
216 1.203, 1.190, 1.173, 1.159, 1.151, 1.145 << 245 {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
217 { 9.714, 7.607, 5.747, 4.493, 3.815, 2.777 << 246 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
218 1.553, 1.353, 1.253, 1.219, 1.211, 1.214 << 247 1.225,1.210,1.191,1.175,1.166,1.174 },
219 1.225, 1.210, 1.191, 1.175, 1.166, 1.174 << 248 {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
220 { 17.97, 12.95, 8.628, 6.065, 4.849, 3.222 << 249 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
221 1.624, 1.382, 1.259, 1.214, 1.202, 1.202 << 250 1.217,1.203,1.184,1.169,1.160,1.151 },
222 1.217, 1.203, 1.184, 1.169, 1.160, 1.151 << 251 {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
223 { 24.83, 17.06, 10.84, 7.355, 5.767, 3.707 << 252 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
224 1.759, 1.465, 1.311, 1.252, 1.234, 1.228 << 253 1.237,1.222,1.201,1.184,1.174,1.159 },
225 1.237, 1.222, 1.201, 1.184, 1.174, 1.159 << 254 {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
226 { 23.26, 17.15, 11.52, 8.049, 6.375, 4.114 << 255 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
227 1.880, 1.535, 1.353, 1.281, 1.258, 1.247 << 256 1.252,1.234,1.212,1.194,1.183,1.170 },
228 1.252, 1.234, 1.212, 1.194, 1.183, 1.170 << 257 {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
229 { 22.33, 18.01, 12.86, 9.212, 7.336, 4.702 << 258 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
230 2.015, 1.602, 1.385, 1.297, 1.268, 1.251 << 259 1.254,1.237,1.214,1.195,1.185,1.179 },
231 1.254, 1.237, 1.214, 1.195, 1.185, 1.179 << 260 {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
232 { 33.91, 24.13, 15.71, 10.80, 8.507, 5.467 << 261 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
233 2.407, 1.873, 1.564, 1.425, 1.374, 1.330 << 262 1.312,1.288,1.258,1.235,1.221,1.205 },
234 1.312, 1.288, 1.258, 1.235, 1.221, 1.205 << 263 {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
235 { 32.14, 24.11, 16.30, 11.40, 9.015, 5.782 << 264 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
236 2.490, 1.925, 1.596, 1.447, 1.391, 1.342 << 265 1.320,1.294,1.264,1.240,1.226,1.214 },
237 1.320, 1.294, 1.264, 1.240, 1.226, 1.214 << 266 {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
238 { 29.51, 24.07, 17.19, 12.28, 9.766, 6.238 << 267 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
239 2.602, 1.995, 1.641, 1.477, 1.414, 1.356 << 268 1.328,1.302,1.270,1.245,1.231,1.233 },
240 1.328, 1.302, 1.270, 1.245, 1.231, 1.233 << 269 {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
241 { 38.19, 30.85, 21.76, 15.35, 12.07, 7.521 << 270 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
242 2.926, 2.188, 1.763, 1.563, 1.484, 1.405 << 271 1.361,1.330,1.294,1.267,1.251,1.239 },
243 1.361, 1.330, 1.294, 1.267, 1.251, 1.239 << 272 {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
244 { 49.71, 39.80, 27.96, 19.63, 15.36, 9.407 << 273 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
245 3.417, 2.478, 1.944, 1.692, 1.589, 1.480 << 274 1.409,1.372,1.330,1.298,1.280,1.258 },
246 1.409, 1.372, 1.330, 1.298, 1.280, 1.258 << 275 {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
247 { 59.25, 45.08, 30.36, 20.83, 16.15, 9.834 << 276 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
248 3.641, 2.648, 2.064, 1.779, 1.661, 1.531 << 277 1.442,1.400,1.354,1.319,1.299,1.272 },
249 1.442, 1.400, 1.354, 1.319, 1.299, 1.272 << 278 {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
250 { 56.38, 44.29, 30.50, 21.18, 16.51, 10.11 << 279 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
251 3.752, 2.724, 2.116, 1.817, 1.692, 1.554 << 280 1.456,1.412,1.364,1.328,1.307,1.282 }};
252 1.456, 1.412, 1.364, 1.328, 1.307, 1.282 << 281
253 }; << 282 //data/corrections for T > Tlim
254 << 283
255 // data/corrections for T > Tlim << 284 static const G4double hecorr[15] = {
256 static constexpr G4double hecorr[15] = { 120 << 285 120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
257 79. << 286 57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
258 41. << 287 -22.30};
259 -7. <<
260 288
261 G4double sigma; 289 G4double sigma;
262 SetParticle(part); 290 SetParticle(part);
263 291
264 Z23 = G4Pow::GetInstance()->Z23(G4lrint(Atom 292 Z23 = G4Pow::GetInstance()->Z23(G4lrint(AtomicNumber));
265 293
266 // correction if particle .ne. e-/e+ 294 // correction if particle .ne. e-/e+
267 // compute equivalent kinetic energy 295 // compute equivalent kinetic energy
268 // lambda depends on p*beta .... 296 // lambda depends on p*beta ....
269 297
270 G4double eKineticEnergy = KineticEnergy; 298 G4double eKineticEnergy = KineticEnergy;
271 299
272 if(mass > electron_mass_c2) 300 if(mass > electron_mass_c2)
273 { 301 {
274 G4double tau1 = KineticEnergy / mass; << 302 G4double TAU = KineticEnergy/mass ;
275 G4double c = mass * tau1 * (tau1 + 2.) / << 303 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
276 G4double w = c - 2.; << 304 G4double w = c-2. ;
277 G4double tau = 0.5 * (w + sqrt(w * w + 4. << 305 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
278 eKineticEnergy = electron_mass_c2 * tau; << 306 eKineticEnergy = electron_mass_c2*tau ;
279 } << 307 }
280 << 308
281 G4double eTotalEnergy = eKineticEnergy + ele << 309 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
282 G4double beta2 = eKineticEnergy * (eT << 310 G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
283 (eTotalEnergy * eTotalEnerg << 311 /(eTotalEnergy*eTotalEnergy);
284 G4double bg2 = eKineticEnergy * (eTotalEnerg << 312 G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
285 (electron_mass_c2 * electron_ << 313 /(electron_mass_c2*electron_mass_c2);
286 << 314
287 static constexpr G4double epsfactor = << 315 static const G4double epsfactor = 2.*CLHEP::electron_mass_c2*
288 2. * CLHEP::electron_mass_c2 * CLHEP::elec << 316 CLHEP::electron_mass_c2*CLHEP::Bohr_radius*CLHEP::Bohr_radius
289 CLHEP::Bohr_radius * CLHEP::Bohr_radius / << 317 /(CLHEP::hbarc*CLHEP::hbarc);
290 G4double eps = epsfactor * bg2 / Z23; << 318 G4double eps = epsfactor*bg2/Z23;
291 << 319
292 if(eps < epsmin) << 320 if (eps<epsmin) sigma = 2.*eps*eps;
293 sigma = 2. * eps * eps; << 321 else if(eps<epsmax) sigma = G4Log(1.+2.*eps)-2.*eps/(1.+2.*eps);
294 else if(eps < epsmax) << 322 else sigma = G4Log(2.*eps)-1.+1./eps;
295 sigma = G4Log(1. + 2. * eps) - 2. * eps / <<
296 else <<
297 sigma = G4Log(2. * eps) - 1. + 1. / eps; <<
298 323
299 sigma *= ChargeSquare * AtomicNumber * Atomi << 324 sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
300 325
301 // interpolate in AtomicNumber and beta2 << 326 // interpolate in AtomicNumber and beta2
302 G4double c1, c2, cc1, cc2, corr; << 327 G4double c1,c2,cc1,cc2,corr;
303 328
304 // get bin number in Z 329 // get bin number in Z
305 G4int iZ = 14; 330 G4int iZ = 14;
306 while((iZ >= 0) && (Zdat[iZ] >= AtomicNumber << 331 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
307 iZ -= 1; << 332 while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
308 if(iZ == 14) << 333 if (iZ==14) iZ = 13;
309 iZ = 13; << 334 if (iZ==-1) iZ = 0 ;
310 if(iZ == -1) <<
311 iZ = 0; <<
312 335
313 G4double ZZ1 = Zdat[iZ]; 336 G4double ZZ1 = Zdat[iZ];
314 G4double ZZ2 = Zdat[iZ + 1]; << 337 G4double ZZ2 = Zdat[iZ+1];
315 G4double ratZ = << 338 G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
316 (AtomicNumber - ZZ1) * (AtomicNumber + ZZ1 << 339 ((ZZ2-ZZ1)*(ZZ2+ZZ1));
317 << 340
318 static constexpr G4double Tlim = 10. * CLHEP << 341 static const G4double Tlim = 10.*CLHEP::MeV;
319 static constexpr G4double sigmafactor = << 342 static const G4double sigmafactor =
320 CLHEP::twopi * CLHEP::classic_electr_radiu << 343 CLHEP::twopi*CLHEP::classic_electr_radius*CLHEP::classic_electr_radius;
321 static const G4double beta2lim = << 344 static const G4double beta2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
322 Tlim * (Tlim + 2. * CLHEP::electron_mass_c << 345 ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHEP::electron_mass_c2));
323 ((Tlim + CLHEP::electron_mass_c2) * (Tlim << 346 static const G4double bg2lim = Tlim*(Tlim+2.*CLHEP::electron_mass_c2)/
324 static const G4double bg2lim = << 347 (CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
325 Tlim * (Tlim + 2. * CLHEP::electron_mass_c << 348
326 (CLHEP::electron_mass_c2 * CLHEP::electron << 349 static const G4double sig0[15] = {
327 << 350 0.2672*CLHEP::barn, 0.5922*CLHEP::barn, 2.653*CLHEP::barn, 6.235*CLHEP::barn,
328 static constexpr G4double sig0[15] = { << 351 11.69*CLHEP::barn , 13.24*CLHEP::barn , 16.12*CLHEP::barn, 23.00*CLHEP::barn,
329 0.2672 * CLHEP::barn, 0.5922 * CLHEP::barn << 352 35.13*CLHEP::barn , 39.95*CLHEP::barn , 50.85*CLHEP::barn, 67.19*CLHEP::barn,
330 6.235 * CLHEP::barn, 11.69 * CLHEP::barn, << 353 91.15*CLHEP::barn , 104.4*CLHEP::barn , 113.1*CLHEP::barn};
331 16.12 * CLHEP::barn, 23.00 * CLHEP::barn, << 354
332 39.95 * CLHEP::barn, 50.85 * CLHEP::barn, << 355 static const G4double Tdat[22] = {
333 91.15 * CLHEP::barn, 104.4 * CLHEP::barn, << 356 100*CLHEP::eV, 200*CLHEP::eV, 400*CLHEP::eV, 700*CLHEP::eV,
334 }; << 357 1*CLHEP::keV, 2*CLHEP::keV, 4*CLHEP::keV, 7*CLHEP::keV,
335 << 358 10*CLHEP::keV, 20*CLHEP::keV, 40*CLHEP::keV, 70*CLHEP::keV,
336 static constexpr G4double Tdat[22] = { << 359 100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP::keV, 700*CLHEP::keV,
337 100 * CLHEP::eV, 200 * CLHEP::eV, 400 * << 360 1*CLHEP::MeV, 2*CLHEP::MeV, 4*CLHEP::MeV, 7*CLHEP::MeV,
338 1 * CLHEP::keV, 2 * CLHEP::keV, 4 * CL << 361 10*CLHEP::MeV, 20*CLHEP::MeV};
339 10 * CLHEP::keV, 20 * CLHEP::keV, 40 * C <<
340 100 * CLHEP::keV, 200 * CLHEP::keV, 400 * <<
341 1 * CLHEP::MeV, 2 * CLHEP::MeV, 4 * CL <<
342 10 * CLHEP::MeV, 20 * CLHEP::MeV <<
343 }; <<
344 362
345 if(eKineticEnergy <= Tlim) << 363 if(eKineticEnergy <= Tlim)
346 { 364 {
347 // get bin number in T (beta2) 365 // get bin number in T (beta2)
348 G4int iT = 21; 366 G4int iT = 21;
349 while((iT >= 0) && (Tdat[iT] >= eKineticEn << 367 // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
350 iT -= 1; << 368 while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
351 if(iT == 21) << 369 if(iT==21) iT = 20;
352 iT = 20; << 370 if(iT==-1) iT = 0 ;
353 if(iT == -1) <<
354 iT = 0; <<
355 371
356 // calculate betasquare values 372 // calculate betasquare values
357 G4double T = Tdat[iT], E = T + electron_ma << 373 G4double T = Tdat[iT], E = T + electron_mass_c2;
358 G4double b2small = T * (E + electron_mass_ << 374 G4double b2small = T*(E+electron_mass_c2)/(E*E);
359 375
360 T = Tdat[iT + 1]; << 376 T = Tdat[iT+1]; E = T + electron_mass_c2;
361 E = T + electron_mass_c2; << 377 G4double b2big = T*(E+electron_mass_c2)/(E*E);
362 G4double b2big = T * (E + electron_mass_c2 << 378 G4double ratb2 = (beta2-b2small)/(b2big-b2small);
363 G4double ratb2 = (beta2 - b2small) / (b2bi <<
364 379
365 if(charge < 0.) << 380 if (charge < 0.)
366 { 381 {
367 c1 = celectron[iZ][iT]; << 382 c1 = celectron[iZ][iT];
368 c2 = celectron[iZ + 1][iT]; << 383 c2 = celectron[iZ+1][iT];
369 cc1 = c1 + ratZ * (c2 - c1); << 384 cc1 = c1+ratZ*(c2-c1);
370 385
371 c1 = celectron[iZ][iT + 1]; << 386 c1 = celectron[iZ][iT+1];
372 c2 = celectron[iZ + 1][iT + 1]; << 387 c2 = celectron[iZ+1][iT+1];
373 cc2 = c1 + ratZ * (c2 - c1); << 388 cc2 = c1+ratZ*(c2-c1);
374 389
375 corr = cc1 + ratb2 * (cc2 - cc1); << 390 corr = cc1+ratb2*(cc2-cc1);
376 391
377 sigma *= sigmafactor / corr; << 392 sigma *= sigmafactor/corr;
378 } 393 }
379 else << 394 else
380 { 395 {
381 c1 = cpositron[iZ][iT]; << 396 c1 = cpositron[iZ][iT];
382 c2 = cpositron[iZ + 1][iT]; << 397 c2 = cpositron[iZ+1][iT];
383 cc1 = c1 + ratZ * (c2 - c1); << 398 cc1 = c1+ratZ*(c2-c1);
384 399
385 c1 = cpositron[iZ][iT + 1]; << 400 c1 = cpositron[iZ][iT+1];
386 c2 = cpositron[iZ + 1][iT + 1]; << 401 c2 = cpositron[iZ+1][iT+1];
387 cc2 = c1 + ratZ * (c2 - c1); << 402 cc2 = c1+ratZ*(c2-c1);
388 403
389 corr = cc1 + ratb2 * (cc2 - cc1); << 404 corr = cc1+ratb2*(cc2-cc1);
390 405
391 sigma *= sigmafactor / corr; << 406 sigma *= sigmafactor/corr;
392 } 407 }
393 } 408 }
394 else 409 else
395 { 410 {
396 c1 = bg2lim * sig0[iZ] * (1. + hecorr[iZ] << 411 c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
397 c2 = << 412 c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
398 bg2lim * sig0[iZ + 1] * (1. + hecorr[iZ <<
399 if((AtomicNumber >= ZZ1) && (AtomicNumber 413 if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
400 sigma = c1 + ratZ * (c2 - c1); << 414 sigma = c1+ratZ*(c2-c1) ;
401 else if(AtomicNumber < ZZ1) 415 else if(AtomicNumber < ZZ1)
402 sigma = AtomicNumber * AtomicNumber * c1 << 416 sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
403 else if(AtomicNumber > ZZ2) 417 else if(AtomicNumber > ZZ2)
404 sigma = AtomicNumber * AtomicNumber * c2 << 418 sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
405 } 419 }
406 return sigma; 420 return sigma;
>> 421
407 } 422 }
408 423
409 //....oooOO0OOooo........oooOO0OOooo........oo 424 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 425
410 void G4UrbanAdjointMscModel::StartTracking(G4T 426 void G4UrbanAdjointMscModel::StartTracking(G4Track* track)
411 { << 427 { SetParticle(track->GetDynamicParticle()->GetDefinition());
412 SetParticle(track->GetDynamicParticle()->Get << 428 firstStep = true;
413 firstStep = true; <<
414 insideskin = false; 429 insideskin = false;
415 fr = facrange; << 430 fr = facrange;
416 tlimit = tgeom = rangeinit = rangecut = geom 431 tlimit = tgeom = rangeinit = rangecut = geombig;
417 smallstep = 1.e1 << 432 smallstep = 1.e10;
418 stepmin = tlim << 433 stepmin = tlimitminfix;
419 tlimitmin = 10. << 434 tlimitmin = 10.*tlimitminfix;
420 rndmEngineMod = G4Ra << 435 rndmEngineMod = G4Random::getTheEngine();
421 } 436 }
422 437
423 //....oooOO0OOooo........oooOO0OOooo........oo 438 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 439
424 G4double G4UrbanAdjointMscModel::ComputeTruePa 440 G4double G4UrbanAdjointMscModel::ComputeTruePathLengthLimit(
425 const G4Track& track, G4double& currentMinim << 441 const G4Track& track,
>> 442 G4double& currentMinimalStep)
426 { 443 {
427 tPathLength = currentMinimal << 444 tPathLength = currentMinimalStep;
428 const G4DynamicParticle* dp = track.GetDynam 445 const G4DynamicParticle* dp = track.GetDynamicParticle();
429 << 446
430 G4StepPoint* sp = track.GetStep()->G << 447 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
431 G4StepStatus stepStatus = sp->GetStepStatus( 448 G4StepStatus stepStatus = sp->GetStepStatus();
432 couple = track.GetMaterialC << 449 couple = track.GetMaterialCutsCouple();
433 SetCurrentCouple(couple); << 450 SetCurrentCouple(couple);
434 currentMaterialIndex = couple->GetIndex(); 451 currentMaterialIndex = couple->GetIndex();
435 currentKinEnergy = dp->GetKineticEnergy( << 452 currentKinEnergy = dp->GetKineticEnergy();
>> 453
>> 454 currentRange = GetRange(particle,currentKinEnergy,couple);
>> 455 lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
>> 456 tPathLength = min(tPathLength,currentRange);
436 457
437 currentRange = GetRange(particle, currentKin <<
438 lambda0 = GetTransportMeanFreePath(part <<
439 tPathLength = min(tPathLength, currentRange <<
440 458
441 // set flag to default values 459 // set flag to default values
442 Zeff = couple->GetMaterial()->GetIonisation( 460 Zeff = couple->GetMaterial()->GetIonisation()->GetZeffective();
>> 461 // couple->GetMaterial()->GetTotNbOfAtomsPerVolume();
443 462
444 if(Zold != Zeff) 463 if(Zold != Zeff)
445 UpdateCache(); 464 UpdateCache();
446 << 465
447 // stop here if small step 466 // stop here if small step
448 if(tPathLength < tlimitminfix) << 467 if(tPathLength < tlimitminfix) {
449 { << 468 latDisplasment = false;
450 latDisplasment = false; << 469 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
451 return ConvertTrueToGeom(tPathLength, curr <<
452 } 470 }
453 471
454 // upper limit for the straight line distanc 472 // upper limit for the straight line distance the particle can travel
455 // for electrons and positrons 473 // for electrons and positrons
456 G4double distance = currentRange; 474 G4double distance = currentRange;
457 // for muons, hadrons 475 // for muons, hadrons
458 if(mass > masslimite) << 476 if(mass > masslimite) {
459 { << 477 distance *= (1.15-9.76e-4*Zeff);
460 distance *= (1.15 - 9.76e-4 * Zeff); << 478 } else {
461 } << 479 distance *= (1.20-Zeff*(1.62e-2-9.22e-5*Zeff));
462 else <<
463 { <<
464 distance *= (1.20 - Zeff * (1.62e-2 - 9.22 <<
465 } 480 }
466 presafety = sp->GetSafety(); 481 presafety = sp->GetSafety();
467 << 482 /*
>> 483 G4cout << "G4Urban::StepLimit tPathLength= "
>> 484 <<tPathLength<<" safety= " << presafety
>> 485 << " range= " <<currentRange<< " lambda= "<<lambda0
>> 486 << " Alg: " << steppingAlgorithm <<G4endl;
>> 487 */
468 // far from geometry boundary 488 // far from geometry boundary
469 if(distance < presafety) 489 if(distance < presafety)
470 { << 490 {
471 latDisplasment = false; << 491 latDisplasment = false;
472 return ConvertTrueToGeom(tPathLength, curr << 492 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
473 } << 493 }
474 494
475 latDisplasment = latDispla << 495 latDisplasment = latDisplasmentbackup;
476 static constexpr G4double invmev = 1.0 / CLH << 496 static const G4double invmev = 1.0/CLHEP::MeV;
>> 497 // standard version
>> 498 //
>> 499 if (steppingAlgorithm == fUseDistanceToBoundary)
>> 500 {
>> 501 //compute geomlimit and presafety
>> 502 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
>> 503 /*
>> 504 G4cout << "G4Urban::Distance to boundary geomlimit= "
>> 505 <<geomlimit<<" safety= " << presafety<<G4endl;
>> 506 */
477 507
478 // standard version << 508 // is it far from boundary ?
479 if(steppingAlgorithm == fUseDistanceToBounda << 509 if(distance < presafety)
480 { << 510 {
481 // compute geomlimit and presafety << 511 latDisplasment = false;
482 geomlimit = ComputeGeomLimit(track, presaf << 512 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 513 }
483 514
484 // is it far from boundary ? << 515 smallstep += 1.;
485 if(distance < presafety) << 516 insideskin = false;
486 { <<
487 latDisplasment = false; <<
488 return ConvertTrueToGeom(tPathLength, cu <<
489 } <<
490 517
491 smallstep += 1.; << 518 // initialisation at firs step and at the boundary
492 insideskin = false; << 519 if(firstStep || (stepStatus == fGeomBoundary))
>> 520 {
>> 521 rangeinit = currentRange;
>> 522 if(!firstStep) { smallstep = 1.; }
493 523
494 // initialisation at firs step and at the << 524 //define stepmin here (it depends on lambda!)
495 if(firstStep || (stepStatus == fGeomBounda << 525 //rough estimation of lambda_elastic/lambda_transport
496 { << 526 G4double rat = currentKinEnergy*invmev;
497 rangeinit = currentRange; << 527 rat = 1.e-3/(rat*(10.+rat)) ;
498 if(!firstStep) << 528 //stepmin ~ lambda_elastic
499 { << 529 stepmin = rat*lambda0;
500 smallstep = 1.; << 530 skindepth = skin*stepmin;
501 } << 531 tlimitmin = max(10*stepmin,tlimitminfix);
>> 532 /*
>> 533 G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
>> 534 << " tlimitmin= " << tlimitmin << " geomlimit= "
>> 535 << geomlimit <<G4endl;
>> 536 */
>> 537 // constraint from the geometry
>> 538
>> 539 if((geomlimit < geombig) && (geomlimit > geommin))
>> 540 {
>> 541 // geomlimit is a geometrical step length
>> 542 // transform it to true path length (estimation)
>> 543 if((1.-geomlimit/lambda0) > 0.)
>> 544 geomlimit = -lambda0*G4Log(1.-geomlimit/lambda0)+tlimitmin ;
>> 545
>> 546 if(stepStatus == fGeomBoundary)
>> 547 tgeom = geomlimit/facgeom;
>> 548 else
>> 549 tgeom = 2.*geomlimit/facgeom;
>> 550 }
>> 551 else
>> 552 tgeom = geombig;
>> 553 }
502 554
503 // define stepmin here (it depends on la << 555 //step limit
504 // rough estimation of lambda_elastic/la << 556 tlimit = facrange*rangeinit;
505 G4double rat = currentKinEnergy * invmev <<
506 rat = 1.e-3 / (rat * (10. + rat <<
507 // stepmin ~ lambda_elastic <<
508 stepmin = rat * lambda0; <<
509 skindepth = skin * stepmin; <<
510 tlimitmin = max(10 * stepmin, tlimitminf <<
511 557
512 // constraint from the geometry << 558 //lower limit for tlimit
513 if((geomlimit < geombig) && (geomlimit > << 559 tlimit = max(tlimit,tlimitmin);
>> 560 tlimit = min(tlimit,tgeom);
>> 561 /*
>> 562 G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
>> 563 << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
>> 564 */
>> 565 // shortcut
>> 566 if((tPathLength < tlimit) && (tPathLength < presafety) &&
>> 567 (smallstep > skin) && (tPathLength < geomlimit-0.999*skindepth))
514 { 568 {
515 // geomlimit is a geometrical step len << 569 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
516 // transform it to true path length (e <<
517 if((1. - geomlimit / lambda0) > 0.) <<
518 geomlimit = -lambda0 * G4Log(1. - ge <<
519 <<
520 if(stepStatus == fGeomBoundary) <<
521 tgeom = geomlimit / facgeom; <<
522 else <<
523 tgeom = 2. * geomlimit / facgeom; <<
524 } 570 }
525 else <<
526 tgeom = geombig; <<
527 } <<
528 571
529 // step limit << 572 // step reduction near to boundary
530 tlimit = facrange * rangeinit; << 573 if(smallstep <= skin)
531 << 574 {
532 // lower limit for tlimit << 575 tlimit = stepmin;
533 tlimit = max(tlimit, tlimitmin); << 576 insideskin = true;
534 tlimit = min(tlimit, tgeom); << 577 }
>> 578 else if(geomlimit < geombig)
>> 579 {
>> 580 if(geomlimit > skindepth)
>> 581 {
>> 582 tlimit = min(tlimit, geomlimit-0.999*skindepth);
>> 583 }
>> 584 else
>> 585 {
>> 586 insideskin = true;
>> 587 tlimit = min(tlimit, stepmin);
>> 588 }
>> 589 }
535 590
536 // shortcut << 591 tlimit = max(tlimit, stepmin);
537 if((tPathLength < tlimit) && (tPathLength <<
538 (smallstep > skin) && (tPathLength < ge <<
539 { <<
540 return ConvertTrueToGeom(tPathLength, cu <<
541 } <<
542 592
543 // step reduction near to boundary << 593 // randomise if not 'small' step and step determined by msc
544 if(smallstep <= skin) << 594 if((tlimit < tPathLength) && (smallstep > skin) && !insideskin)
545 { << 595 {
546 tlimit = stepmin; << 596 tPathLength = min(tPathLength, Randomizetlimit());
547 insideskin = true; << 597 }
548 } <<
549 else if(geomlimit < geombig) <<
550 { <<
551 if(geomlimit > skindepth) <<
552 { <<
553 tlimit = min(tlimit, geomlimit - 0.999 <<
554 } <<
555 else 598 else
556 { << 599 {
557 insideskin = true; << 600 tPathLength = min(tPathLength, tlimit);
558 tlimit = min(tlimit, stepmin); << 601 }
559 } <<
560 } <<
561 <<
562 tlimit = max(tlimit, stepmin); <<
563 602
564 // randomise if not 'small' step and step <<
565 if((tlimit < tPathLength) && (smallstep > <<
566 { <<
567 tPathLength = min(tPathLength, Randomize <<
568 } 603 }
569 else << 604 // for 'normal' simulation with or without magnetic field
570 { << 605 // there no small step/single scattering at boundaries
571 tPathLength = min(tPathLength, tlimit); <<
572 } <<
573 } <<
574 // for 'normal' simulation with or without m <<
575 // there no small step/single scattering at <<
576 else if(steppingAlgorithm == fUseSafety) 606 else if(steppingAlgorithm == fUseSafety)
577 { <<
578 if(stepStatus != fGeomBoundary) <<
579 { <<
580 presafety = ComputeSafety(sp->GetPositio <<
581 } <<
582 // is far from boundary <<
583 if(distance < presafety) <<
584 { 607 {
585 latDisplasment = false; << 608 if(stepStatus != fGeomBoundary) {
586 return ConvertTrueToGeom(tPathLength, cu << 609 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
587 } << 610 }
588 << 611 /*
589 if(firstStep || (stepStatus == fGeomBounda << 612 G4cout << "presafety= " << presafety
590 { << 613 << " firstStep= " << firstStep
591 rangeinit = currentRange; << 614 << " stepStatus= " << stepStatus
592 fr = facrange; << 615 << G4endl;
593 // Geant4 version 9.1 like stepping for << 616 */
594 if(mass < masslimite) << 617 // is far from boundary
595 { << 618 if(distance < presafety)
596 rangeinit = max(rangeinit, lambda0); <<
597 if(lambda0 > lambdalimit) <<
598 { 619 {
599 fr *= (0.75 + 0.25 * lambda0 / lambd << 620 latDisplasment = false;
>> 621 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
600 } 622 }
601 } <<
602 // lower limit for tlimit <<
603 G4double rat = currentKinEnergy * invmev <<
604 rat = 1.e-3 / (rat * (10 + rat) <<
605 stepmin = lambda0 * rat; <<
606 tlimitmin = max(10 * stepmin, tlimitm <<
607 } <<
608 <<
609 // step limit <<
610 tlimit = max(fr * rangeinit, facsafety * p <<
611 623
612 // lower limit for tlimit << 624 if(firstStep || (stepStatus == fGeomBoundary)) {
613 tlimit = max(tlimit, tlimitmin); << 625 rangeinit = currentRange;
614 << 626 fr = facrange;
615 // randomise if step determined by msc << 627 // 9.1 like stepping for e+/e- only (not for muons,hadrons)
616 if(tlimit < tPathLength) << 628 if(mass < masslimite)
617 { << 629 {
618 tPathLength = min(tPathLength, Randomize << 630 rangeinit = max(rangeinit, lambda0);
619 } << 631 if(lambda0 > lambdalimit) {
620 else << 632 fr *= (0.75+0.25*lambda0/lambdalimit);
621 { << 633 }
622 tPathLength = min(tPathLength, tlimit); << 634 }
>> 635 //lower limit for tlimit
>> 636 G4double rat = currentKinEnergy*invmev;
>> 637 rat = 1.e-3/(rat*(10 + rat)) ;
>> 638 stepmin = lambda0*rat;
>> 639 tlimitmin = max(10*stepmin, tlimitminfix);
>> 640 }
>> 641
>> 642 //step limit
>> 643 tlimit = max(fr*rangeinit, facsafety*presafety);
>> 644
>> 645 //lower limit for tlimit
>> 646 tlimit = max(tlimit, tlimitmin);
>> 647
>> 648 // randomise if step determined by msc
>> 649 if(tlimit < tPathLength)
>> 650 {
>> 651 tPathLength = min(tPathLength, Randomizetlimit());
>> 652 }
>> 653 else { tPathLength = min(tPathLength, tlimit); }
623 } 654 }
624 } <<
625 // new stepping mode UseSafetyPlus 655 // new stepping mode UseSafetyPlus
626 else if(steppingAlgorithm == fUseSafetyPlus) 656 else if(steppingAlgorithm == fUseSafetyPlus)
627 { <<
628 if(stepStatus != fGeomBoundary) <<
629 { 657 {
630 presafety = ComputeSafety(sp->GetPositio << 658 if(stepStatus != fGeomBoundary) {
631 } << 659 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
632 << 660 }
633 // is far from boundary << 661 /*
634 if(distance < presafety) << 662 G4cout << "presafety= " << presafety
635 { << 663 << " firstStep= " << firstStep
636 latDisplasment = false; << 664 << " stepStatus= " << stepStatus
637 return ConvertTrueToGeom(tPathLength, cu << 665 << G4endl;
638 } << 666 */
639 << 667 // is far from boundary
640 if(firstStep || (stepStatus == fGeomBounda << 668 if(distance < presafety)
641 { <<
642 rangeinit = currentRange; <<
643 fr = facrange; <<
644 rangecut = geombig; <<
645 if(mass < masslimite) <<
646 { <<
647 G4int index = 1; <<
648 if(charge > 0.) <<
649 index = 2; <<
650 rangecut = couple->GetProductionCuts() <<
651 if(lambda0 > lambdalimit) <<
652 { 669 {
653 fr *= (0.84 + 0.16 * lambda0 / lambd << 670 latDisplasment = false;
>> 671 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
654 } 672 }
>> 673
>> 674 if(firstStep || (stepStatus == fGeomBoundary)) {
>> 675 rangeinit = currentRange;
>> 676 fr = facrange;
>> 677 rangecut = geombig;
>> 678 if(mass < masslimite)
>> 679 {
>> 680 G4int index = 1;
>> 681 if(charge > 0.) index = 2;
>> 682 rangecut = couple->GetProductionCuts()->GetProductionCut(index);
>> 683 if(lambda0 > lambdalimit) {
>> 684 fr *= (0.84+0.16*lambda0/lambdalimit);
>> 685 }
>> 686 }
>> 687 //lower limit for tlimit
>> 688 G4double rat = currentKinEnergy*invmev;
>> 689 rat = 1.e-3/(rat*(10 + rat)) ;
>> 690 stepmin = lambda0*rat;
>> 691 tlimitmin = max(10*stepmin, tlimitminfix);
655 } 692 }
656 // lower limit for tlimit << 693 //step limit
657 G4double rat = currentKinEnergy * invmev << 694 tlimit = max(fr*rangeinit, facsafety*presafety);
658 rat = 1.e-3 / (rat * (10 + rat) <<
659 stepmin = lambda0 * rat; <<
660 tlimitmin = max(10 * stepmin, tlimitm <<
661 } <<
662 // step limit <<
663 tlimit = max(fr * rangeinit, facsafety * p <<
664 695
665 // lower limit for tlimit << 696 //lower limit for tlimit
666 tlimit = max(tlimit, tlimitmin); << 697 tlimit = max(tlimit, tlimitmin);
667 698
668 // condition for tPathLength from drr and << 699 // condition for tPathLength from drr and finalr
669 if(currentRange > finalr) << 700 if(currentRange > finalr) {
670 { << 701 G4double tmax = drr*currentRange+
671 G4double tmax = << 702 finalr*(1.-drr)*(2.-finalr/currentRange);
672 drr * currentRange + finalr * (1. - dr << 703 tPathLength = min(tPathLength,tmax);
673 tPathLength = min(tPathLength, tmax); << 704 }
674 } <<
675 705
676 // condition safety << 706 // condition safety
677 if(currentRange > rangecut) << 707 if(currentRange > rangecut) {
678 { << 708 if(firstStep) {
679 if(firstStep) << 709 tPathLength = min(tPathLength,facsafety*presafety);
680 { << 710 } else if(stepStatus != fGeomBoundary && presafety > stepmin) {
681 tPathLength = min(tPathLength, facsafe << 711 tPathLength = min(tPathLength,presafety);
>> 712 }
682 } 713 }
683 else if(stepStatus != fGeomBoundary && p << 714
>> 715 // randomise if step determined by msc
>> 716 if(tPathLength < tlimit)
684 { 717 {
685 tPathLength = min(tPathLength, presafe << 718 tPathLength = min(tPathLength, Randomizetlimit());
686 } 719 }
>> 720 else { tPathLength = min(tPathLength, tlimit); }
687 } 721 }
688 722
689 // randomise if step determined by msc <<
690 if(tPathLength < tlimit) <<
691 { <<
692 tPathLength = min(tPathLength, Randomize <<
693 } <<
694 else <<
695 { <<
696 tPathLength = min(tPathLength, tlimit); <<
697 } <<
698 } <<
699 <<
700 // version similar to 7.1 (needed for some e 723 // version similar to 7.1 (needed for some experiments)
701 else 724 else
702 { <<
703 if(stepStatus == fGeomBoundary) <<
704 { 725 {
705 if(currentRange > lambda0) << 726 if (stepStatus == fGeomBoundary)
706 { << 727 {
707 tlimit = facrange * currentRange; << 728 if (currentRange > lambda0) { tlimit = facrange*currentRange; }
708 } << 729 else { tlimit = facrange*lambda0; }
709 else << 730
>> 731 tlimit = max(tlimit, tlimitmin);
>> 732 }
>> 733 // randomise if step determined by msc
>> 734 if(tlimit < tPathLength)
710 { 735 {
711 tlimit = facrange * lambda0; << 736 tPathLength = min(tPathLength, Randomizetlimit());
712 } 737 }
713 << 738 else { tPathLength = min(tPathLength, tlimit); }
714 tlimit = max(tlimit, tlimitmin); <<
715 } 739 }
716 // randomise if step determined by msc << 740 firstStep = false;
717 if(tlimit < tPathLength) <<
718 { <<
719 tPathLength = min(tPathLength, Randomize <<
720 } <<
721 else <<
722 { <<
723 tPathLength = min(tPathLength, tlimit); <<
724 } <<
725 } <<
726 firstStep = false; <<
727 return ConvertTrueToGeom(tPathLength, curren 741 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
728 } 742 }
729 743
730 //....oooOO0OOooo........oooOO0OOooo........oo 744 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 745
731 G4double G4UrbanAdjointMscModel::ComputeGeomPa 746 G4double G4UrbanAdjointMscModel::ComputeGeomPathLength(G4double)
732 { 747 {
733 lambdaeff = lambda0; 748 lambdaeff = lambda0;
734 par1 = -1.; << 749 par1 = -1. ;
735 par2 = par3 = 0.; << 750 par2 = par3 = 0. ;
736 751
737 // this correction needed to run MSC with eI 752 // this correction needed to run MSC with eIoni and eBrem inactivated
738 // and makes no harm for a normal run 753 // and makes no harm for a normal run
739 tPathLength = std::min(tPathLength, currentR << 754 tPathLength = std::min(tPathLength,currentRange);
740 755
741 // do the true -> geom transformation 756 // do the true -> geom transformation
742 zPathLength = tPathLength; 757 zPathLength = tPathLength;
743 758
744 // z = t for very small tPathLength 759 // z = t for very small tPathLength
745 if(tPathLength < tlimitminfix2) << 760 if(tPathLength < tlimitminfix2) return zPathLength;
746 return zPathLength; <<
747 <<
748 G4double tau = tPathLength / lambda0; <<
749 761
750 if((tau <= tausmall) || insideskin) << 762 // VI: it is already checked
751 { << 763 // if(tPathLength > currentRange)
752 zPathLength = min(tPathLength, lambda0); << 764 // tPathLength = currentRange ;
753 } << 765 /*
754 else if(tPathLength < currentRange * dtrl) << 766 G4cout << "ComputeGeomPathLength: tpl= " << tPathLength
755 { << 767 << " R= " << currentRange << " L0= " << lambda0
756 if(tau < taulim) << 768 << " E= " << currentKinEnergy << " "
757 zPathLength = tPathLength * (1. - 0.5 * << 769 << particle->GetParticleName() << G4endl;
758 else << 770 */
759 zPathLength = lambda0 * (1. - G4Exp(-tau << 771 G4double tau = tPathLength/lambda0 ;
760 } << 772
761 else if(currentKinEnergy < mass || tPathLeng << 773 if ((tau <= tausmall) || insideskin) {
762 { << 774 zPathLength = min(tPathLength, lambda0);
763 par1 = 1. / currentRange; << 775
764 par2 = 1. / (par1 * lambda0); << 776 } else if (tPathLength < currentRange*dtrl) {
765 par3 = 1. + par2; << 777 if(tau < taulim) zPathLength = tPathLength*(1.-0.5*tau) ;
766 if(tPathLength < currentRange) << 778 else zPathLength = lambda0*(1.-G4Exp(-tau));
767 { << 779
768 zPathLength = << 780 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
769 (1. - G4Exp(par3 * G4Log(1. - tPathLen << 781 par1 = 1./currentRange ;
770 (par1 * par3); << 782 par2 = 1./(par1*lambda0) ;
771 } << 783 par3 = 1.+par2 ;
772 else << 784 if(tPathLength < currentRange) {
773 { << 785 zPathLength =
774 zPathLength = 1. / (par1 * par3); << 786 (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3);
775 } << 787 } else {
776 } << 788 zPathLength = 1./(par1*par3);
777 else << 789 }
778 { << 790
779 G4double rfin = max(currentRange - tPat << 791 } else {
780 G4double T1 = GetEnergy(particle, rfi << 792 G4double rfin = max(currentRange-tPathLength, 0.01*currentRange);
781 G4double lambda1 = GetTransportMeanFreePat << 793 G4double T1 = GetEnergy(particle,rfin,couple);
782 << 794 G4double lambda1 = GetTransportMeanFreePath(particle,T1);
783 par1 = (lambda0 - lambda1) / (lambd << 795
784 par2 = 1. / (par1 * lambda0); << 796 par1 = (lambda0-lambda1)/(lambda0*tPathLength);
785 par3 = 1. + par2; << 797 //G4cout << "par1= " << par1 << " L1= " << lambda1 << G4endl;
786 zPathLength = (1. - G4Exp(par3 * G4Log(lam << 798 par2 = 1./(par1*lambda0);
>> 799 par3 = 1.+par2 ;
>> 800 zPathLength = (1.-G4Exp(par3*G4Log(lambda1/lambda0)))/(par1*par3);
787 } 801 }
788 802
789 zPathLength = min(zPathLength, lambda0); 803 zPathLength = min(zPathLength, lambda0);
>> 804 //G4cout<< "zPathLength= "<< zPathLength<< " L0= " << lambda0 << G4endl;
790 return zPathLength; 805 return zPathLength;
791 } 806 }
792 807
793 //....oooOO0OOooo........oooOO0OOooo........oo 808 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 809
794 G4double G4UrbanAdjointMscModel::ComputeTrueSt 810 G4double G4UrbanAdjointMscModel::ComputeTrueStepLength(G4double geomStepLength)
795 { 811 {
796 // step defined other than transportation << 812 // step defined other than transportation
797 if(geomStepLength == zPathLength) << 813 if(geomStepLength == zPathLength) {
798 { << 814 //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
799 return tPathLength; << 815 // << " step= " << geomStepLength << " *** " << G4endl;
>> 816 return tPathLength;
800 } 817 }
801 818
802 zPathLength = geomStepLength; 819 zPathLength = geomStepLength;
803 820
804 // t = z for very small step 821 // t = z for very small step
805 if(geomStepLength < tlimitminfix2) << 822 if(geomStepLength < tlimitminfix2) {
806 { << 823 tPathLength = geomStepLength;
807 tPathLength = geomStepLength; << 824
>> 825 // recalculation
>> 826 } else {
808 827
809 // recalculation <<
810 } <<
811 else <<
812 { <<
813 G4double tlength = geomStepLength; 828 G4double tlength = geomStepLength;
814 if((geomStepLength > lambda0 * tausmall) & << 829 if((geomStepLength > lambda0*tausmall) && !insideskin) {
815 { << 830
816 if(par1 < 0.) << 831 if(par1 < 0.) {
817 { << 832 tlength = -lambda0*G4Log(1.-geomStepLength/lambda0) ;
818 tlength = -lambda0 * G4Log(1. - geomSt << 833 } else {
819 } << 834 if(par1*par3*geomStepLength < 1.) {
820 else << 835 tlength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
821 { << 836 } else {
822 if(par1 * par3 * geomStepLength < 1.) <<
823 { <<
824 tlength = <<
825 (1. - G4Exp(G4Log(1. - par1 * par3 <<
826 par1; <<
827 } <<
828 else <<
829 { <<
830 tlength = currentRange; 837 tlength = currentRange;
831 } 838 }
832 } 839 }
833 840
834 if(tlength < geomStepLength) << 841 if(tlength < geomStepLength) { tlength = geomStepLength; }
835 { << 842 else if(tlength > tPathLength) { tlength = tPathLength; }
836 tlength = geomStepLength; << 843 }
837 } << 844 tPathLength = tlength;
838 else if(tlength > tPathLength) <<
839 { <<
840 tlength = tPathLength; <<
841 } <<
842 } <<
843 tPathLength = tlength; <<
844 } 845 }
>> 846 //G4cout << "Urban::ComputeTrueLength: tPathLength= " << tPathLength
>> 847 // << " step= " << geomStepLength << " &&& " << G4endl;
845 848
846 return tPathLength; 849 return tPathLength;
847 } 850 }
848 851
849 //....oooOO0OOooo........oooOO0OOooo........oo 852 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
850 G4ThreeVector& G4UrbanAdjointMscModel::SampleS << 853
851 const G4ThreeVector& oldDirection, G4double << 854 G4ThreeVector&
>> 855 G4UrbanAdjointMscModel::SampleScattering(const G4ThreeVector& oldDirection,
>> 856 G4double /*safety*/)
852 { 857 {
853 fDisplacement.set(0.0, 0.0, 0.0); << 858 fDisplacement.set(0.0,0.0,0.0);
854 G4double kineticEnergy = currentKinEnergy; 859 G4double kineticEnergy = currentKinEnergy;
855 if(tPathLength > currentRange * dtrl) << 860 if (tPathLength > currentRange*dtrl) {
856 { << 861 kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
857 kineticEnergy = GetEnergy(particle, curren << 862 } else {
858 } << 863 kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
859 else <<
860 { <<
861 kineticEnergy -= tPathLength * GetDEDX(par <<
862 } 864 }
863 865
864 if((kineticEnergy <= eV) || (tPathLength <= 866 if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
865 (tPathLength < tausmall * lambda0)) << 867 (tPathLength < tausmall*lambda0)) { return fDisplacement; }
866 { <<
867 return fDisplacement; <<
868 } <<
869 868
870 G4double cth = SampleCosineTheta(tPathLength << 869 G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
871 870
872 // protection against 'bad' cth values 871 // protection against 'bad' cth values
873 if(std::fabs(cth) >= 1.0) << 872 if(std::fabs(cth) >= 1.0) { return fDisplacement; }
874 { <<
875 return fDisplacement; <<
876 } <<
877 873
878 G4double sth = sqrt((1.0 - cth) * (1.0 + ct << 874 /*
879 G4double phi = twopi * rndmEngineMod->flat( << 875 if(cth < 1.0 - 1000*tPathLength/lambda0 && cth < 0.5 &&
880 G4double dirx = sth * cos(phi); << 876 kineticEnergy > 20*MeV) {
881 G4double diry = sth * sin(phi); << 877 G4cout << "### G4UrbanAdjointMscModel::SampleScattering for "
>> 878 << particle->GetParticleName()
>> 879 << " E(MeV)= " << kineticEnergy/MeV
>> 880 << " Step(mm)= " << tPathLength/mm
>> 881 << " in " << CurrentCouple()->GetMaterial()->GetName()
>> 882 << " CosTheta= " << cth << G4endl;
>> 883 }
>> 884 */
>> 885 G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
>> 886 G4double phi = twopi*rndmEngineMod->flat();
>> 887 G4double dirx = sth*cos(phi);
>> 888 G4double diry = sth*sin(phi);
882 889
883 G4ThreeVector newDirection(dirx, diry, cth); << 890 G4ThreeVector newDirection(dirx,diry,cth);
884 newDirection.rotateUz(oldDirection); 891 newDirection.rotateUz(oldDirection);
885 892
886 fParticleChange->ProposeMomentumDirection(ne 893 fParticleChange->ProposeMomentumDirection(newDirection);
>> 894 /*
>> 895 G4cout << "G4UrbanAdjointMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
>> 896 << " sinTheta= " << sth << " safety(mm)= " << safety
>> 897 << " trueStep(mm)= " << tPathLength
>> 898 << " geomStep(mm)= " << zPathLength
>> 899 << G4endl;
>> 900 */
887 901
888 if(latDisplasment && currentTau >= tausmall) << 902
889 { << 903 if (latDisplasment && currentTau >= tausmall) {
890 if(displacementFlag) << 904 if(displacementFlag) { SampleDisplacementNew(cth, phi); }
891 { << 905 else { SampleDisplacement(sth, phi); }
892 SampleDisplacementNew(cth, phi); <<
893 } <<
894 else <<
895 { <<
896 SampleDisplacement(sth, phi); <<
897 } <<
898 fDisplacement.rotateUz(oldDirection); 906 fDisplacement.rotateUz(oldDirection);
899 } 907 }
900 return fDisplacement; 908 return fDisplacement;
901 } 909 }
902 910
903 //....oooOO0OOooo........oooOO0OOooo........oo 911 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 912
904 G4double G4UrbanAdjointMscModel::SampleCosineT 913 G4double G4UrbanAdjointMscModel::SampleCosineTheta(G4double trueStepLength,
905 << 914 G4double KineticEnergy)
906 { 915 {
907 G4double cth = 1.; << 916 G4double cth = 1. ;
908 G4double tau = trueStepLength / lambda0; << 917 G4double tau = trueStepLength/lambda0;
909 currentTau = tau; 918 currentTau = tau;
910 lambdaeff = lambda0; 919 lambdaeff = lambda0;
911 920
912 G4double lambda1 = GetTransportMeanFreePath( << 921 G4double lambda1 = GetTransportMeanFreePath(particle,KineticEnergy);
913 if(std::fabs(lambda1 - lambda0) > lambda0 * << 922 if(std::fabs(lambda1 - lambda0) > lambda0*0.01 && lambda1 > 0.)
914 { 923 {
915 // mean tau value 924 // mean tau value
916 tau = trueStepLength * G4Log(lambda0 / lam << 925 tau = trueStepLength*G4Log(lambda0/lambda1)/(lambda0-lambda1);
917 } 926 }
918 927
919 currentTau = tau; << 928 currentTau = tau ;
920 lambdaeff = trueStepLength / currentT << 929 lambdaeff = trueStepLength/currentTau;
921 currentRadLength = couple->GetMaterial()->Ge 930 currentRadLength = couple->GetMaterial()->GetRadlen();
922 931
923 if(tau >= taubig) << 932 if (tau >= taubig) { cth = -1.+2.*rndmEngineMod->flat(); }
924 { << 933 else if (tau >= tausmall) {
925 cth = -1. + 2. * rndmEngineMod->flat(); <<
926 } <<
927 else if(tau >= tausmall) <<
928 { <<
929 static const G4double numlim = 0.01; 934 static const G4double numlim = 0.01;
930 G4double xmeanth, x2meanth; 935 G4double xmeanth, x2meanth;
931 if(tau < numlim) << 936 if(tau < numlim) {
932 { << 937 xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
933 xmeanth = 1.0 - tau * (1.0 - 0.5 * tau) << 938 x2meanth= 1.0 - tau*(5.0 - 6.25*tau)/3.;
934 x2meanth = 1.0 - tau * (5.0 - 6.25 * tau << 939 } else {
935 } << 940 xmeanth = G4Exp(-tau);
936 else << 941 x2meanth = (1.+2.*G4Exp(-2.5*tau))/3.;
937 { <<
938 xmeanth = G4Exp(-tau); <<
939 x2meanth = (1. + 2. * G4Exp(-2.5 * tau)) <<
940 } 942 }
941 943
942 // too large step of low-energy particle 944 // too large step of low-energy particle
943 G4double relloss = 1. - Ki << 945 G4double relloss = 1. - KineticEnergy/currentKinEnergy;
944 static const G4double rellossmax = 0.50; << 946 static const G4double rellossmax= 0.50;
945 if(relloss > rellossmax) << 947 if(relloss > rellossmax) {
946 { << 948 return SimpleScattering(xmeanth,x2meanth);
947 return SimpleScattering(xmeanth, x2meant <<
948 } 949 }
949 // is step extreme small ? 950 // is step extreme small ?
950 G4bool extremesmallstep = false; << 951 G4bool extremesmallstep = false ;
951 G4double tsmall = std::min(tlimitm << 952 G4double tsmall = std::min(tlimitmin,lambdalimit);
952 G4double theta0 = 0.; << 953 G4double theta0 = 0.;
953 if(trueStepLength > tsmall) << 954 if(trueStepLength > tsmall) {
954 { << 955 theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
955 theta0 = ComputeTheta0(trueStepLength, K << 956 } else {
956 } << 957 theta0 = sqrt(trueStepLength/tsmall)*ComputeTheta0(tsmall,KineticEnergy);
957 else << 958 extremesmallstep = true ;
958 { <<
959 theta0 = <<
960 sqrt(trueStepLength / tsmall) * Comput <<
961 extremesmallstep = true; <<
962 } 959 }
963 960
964 static constexpr G4double theta0max = CLHE << 961 static const G4double theta0max = CLHEP::pi/6.;
>> 962 //G4cout << "Theta0= " << theta0 << " theta0max= " << theta0max
>> 963 // << " sqrt(tausmall)= " << sqrt(tausmall) << G4endl;
965 964
966 // protection for very small angles 965 // protection for very small angles
967 G4double theta2 = theta0 * theta0; << 966 G4double theta2 = theta0*theta0;
968 967
969 if(theta2 < tausmall) << 968 if(theta2 < tausmall) { return cth; }
970 { << 969
971 return cth; << 970 if(theta0 > theta0max) {
>> 971 return SimpleScattering(xmeanth,x2meanth);
972 } 972 }
973 973
974 if(theta0 > theta0max) << 974 G4double x = theta2*(1.0 - theta2/12.);
975 { << 975 if(theta2 > numlim) {
976 return SimpleScattering(xmeanth, x2meant << 976 G4double sth = 2*sin(0.5*theta0);
977 } << 977 x = sth*sth;
978 <<
979 G4double x = theta2 * (1.0 - theta2 / 12.) <<
980 if(theta2 > numlim) <<
981 { <<
982 G4double sth = 2 * sin(0.5 * theta0); <<
983 x = sth * sth; <<
984 } 978 }
985 979
986 // parameter for tail 980 // parameter for tail
987 G4double ltau = G4Log(tau); << 981 G4double ltau= G4Log(tau);
988 G4double u = G4Exp(ltau / 6.); << 982 G4double u = G4Exp(ltau/6.);
989 if(extremesmallstep) << 983 if(extremesmallstep) u = G4Exp(G4Log(tsmall/lambda0)/6.);
990 u = G4Exp(G4Log(tsmall / lambda0) / 6.); << 984 G4double xx = G4Log(lambdaeff/currentRadLength);
991 G4double xx = G4Log(lambdaeff / currentRa << 985 G4double xsi = coeffc1+u*(coeffc2+coeffc3*u)+coeffc4*xx;
992 G4double xsi = coeffc1 + u * (coeffc2 + co <<
993 986
994 // tail should not be too big 987 // tail should not be too big
995 if(xsi < 1.9) << 988 if(xsi < 1.9) {
996 { << 989 /*
997 xsi = 1.9; << 990 if(KineticEnergy > 20*MeV && xsi < 1.6) {
>> 991 G4cout << "G4UrbanAdjointMscModel::SampleCosineTheta: E(GeV)= "
>> 992 << KineticEnergy/GeV
>> 993 << " !!** c= " << xsi
>> 994 << " **!! length(mm)= " << trueStepLength << " Zeff= " << Zeff
>> 995 << " " << couple->GetMaterial()->GetName()
>> 996 << " tau= " << tau << G4endl;
>> 997 }
>> 998 */
>> 999 xsi = 1.9;
998 } 1000 }
999 1001
1000 G4double c = xsi; 1002 G4double c = xsi;
1001 1003
1002 if(std::abs(c - 3.) < 0.001) << 1004 if(std::abs(c-3.) < 0.001) { c = 3.001; }
1003 { << 1005 else if(std::abs(c-2.) < 0.001) { c = 2.001; }
1004 c = 3.001; <<
1005 } <<
1006 else if(std::abs(c - 2.) < 0.001) <<
1007 { <<
1008 c = 2.001; <<
1009 } <<
1010 1006
1011 G4double c1 = c - 1.; << 1007 G4double c1 = c-1.;
1012 1008
1013 G4double ea = G4Exp(-xsi); << 1009 G4double ea = G4Exp(-xsi);
1014 G4double eaa = 1. - ea; << 1010 G4double eaa = 1.-ea ;
1015 G4double xmean1 = 1. - (1. - (1. + xsi) * << 1011 G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
1016 G4double x0 = 1. - xsi * x; << 1012 G4double x0 = 1. - xsi*x;
1017 1013
1018 if(xmean1 <= 0.999 * xmeanth) << 1014 // G4cout << " xmean1= " << xmean1 << " xmeanth= " << xmeanth << G4endl;
1019 { << 1015
1020 return SimpleScattering(xmeanth, x2mean << 1016 if(xmean1 <= 0.999*xmeanth) {
>> 1017 return SimpleScattering(xmeanth,x2meanth);
1021 } 1018 }
1022 // from continuity of derivatives << 1019 //from continuity of derivatives
1023 G4double b = 1. + (c - xsi) * x; << 1020 G4double b = 1.+(c-xsi)*x;
1024 1021
1025 G4double b1 = b + 1.; << 1022 G4double b1 = b+1.;
1026 G4double bx = c * x; << 1023 G4double bx = c*x;
1027 1024
1028 G4double eb1 = G4Exp(G4Log(b1) * c1); << 1025 G4double eb1 = G4Exp(G4Log(b1)*c1);
1029 G4double ebx = G4Exp(G4Log(bx) * c1); << 1026 G4double ebx = G4Exp(G4Log(bx)*c1);
1030 G4double d = ebx / eb1; << 1027 G4double d = ebx/eb1;
1031 1028
1032 G4double xmean2 = (x0 + d - (bx - b1 * d) << 1029 G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
1033 1030
1034 G4double f1x0 = ea / eaa; << 1031 G4double f1x0 = ea/eaa;
1035 G4double f2x0 = c1 / (c * (1. - d)); << 1032 G4double f2x0 = c1/(c*(1. - d));
1036 G4double prob = f2x0 / (f1x0 + f2x0); << 1033 G4double prob = f2x0/(f1x0+f2x0);
1037 1034
1038 G4double qprob = xmeanth / (prob * xmean1 << 1035 G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1039 1036
>> 1037 // sampling of costheta
>> 1038 //G4cout << "c= " << c << " qprob= " << qprob << " eb1= " << eb1
>> 1039 // << " c1= " << c1 << " b1= " << b1 << " bx= " << bx << " eb1= " << eb1
>> 1040 // << G4endl;
1040 if(rndmEngineMod->flat() < qprob) 1041 if(rndmEngineMod->flat() < qprob)
1041 { 1042 {
1042 G4double var = 0; 1043 G4double var = 0;
1043 if(rndmEngineMod->flat() < prob) << 1044 if(rndmEngineMod->flat() < prob) {
1044 { << 1045 cth = 1.+G4Log(ea+rndmEngineMod->flat()*eaa)*x;
1045 cth = 1. + G4Log(ea + rndmEngineMod-> << 1046 } else {
1046 } << 1047 var = (1.0 - d)*rndmEngineMod->flat();
1047 else << 1048 if(var < numlim*d) {
1048 { << 1049 var /= (d*c1);
1049 var = (1.0 - d) * rndmEngineMod->flat << 1050 cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
1050 if(var < numlim * d) << 1051 } else {
1051 { << 1052 cth = 1. + x*(c - xsi - c*G4Exp(-G4Log(var + d)/c1));
1052 var /= (d * c1); <<
1053 cth = -1.0 + var * (1.0 - 0.5 * var <<
1054 } <<
1055 else <<
1056 { <<
1057 cth = 1. + x * (c - xsi - c * G4Exp <<
1058 } 1053 }
>> 1054 }
>> 1055 /*
>> 1056 if(KineticEnergy > 5*GeV && cth < 0.9) {
>> 1057 G4cout << "G4UrbanAdjointMscModel::SampleCosineTheta: E(GeV)= "
>> 1058 << KineticEnergy/GeV
>> 1059 << " 1-cosT= " << 1 - cth
>> 1060 << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
>> 1061 << " tau= " << tau
>> 1062 << " prob= " << prob << " var= " << var << G4endl;
>> 1063 G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
>> 1064 << " ebx= " << ebx
>> 1065 << " c1= " << c1 << " b= " << b << " b1= " << b1
>> 1066 << " bx= " << bx << " d= " << d
>> 1067 << " ea= " << ea << " eaa= " << eaa << G4endl;
>> 1068 }
>> 1069 */
>> 1070 }
>> 1071 else {
>> 1072 cth = -1.+2.*rndmEngineMod->flat();
>> 1073 /*
>> 1074 if(KineticEnergy > 5*GeV) {
>> 1075 G4cout << "G4UrbanAdjointMscModel::SampleCosineTheta: E(GeV)= "
>> 1076 << KineticEnergy/GeV
>> 1077 << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
>> 1078 << " qprob= " << qprob << G4endl;
1059 } 1079 }
1060 } << 1080 */
1061 else <<
1062 { <<
1063 cth = -1. + 2. * rndmEngineMod->flat(); <<
1064 } 1081 }
1065 } 1082 }
1066 return cth; << 1083 return cth ;
1067 } 1084 }
1068 1085
1069 //....oooOO0OOooo........oooOO0OOooo........o 1086 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1087
1070 G4double G4UrbanAdjointMscModel::ComputeTheta 1088 G4double G4UrbanAdjointMscModel::ComputeTheta0(G4double trueStepLength,
1071 << 1089 G4double KineticEnergy)
1072 { 1090 {
1073 // for all particles take the width of the 1091 // for all particles take the width of the central part
1074 // from a parametrization similar to the H << 1092 // from a parametrization similar to the Highland formula
1075 // ( Highland formula: Particle Physics Boo 1093 // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
1076 G4double invbetacp = << 1094 G4double invbetacp = std::sqrt((currentKinEnergy+mass)*(KineticEnergy+mass)/
1077 std::sqrt((currentKinEnergy + mass) * (Ki << 1095 (currentKinEnergy*(currentKinEnergy+2.*mass)*
1078 (currentKinEnergy * (currentKin << 1096 KineticEnergy*(KineticEnergy+2.*mass)));
1079 KineticEnergy * (KineticEnergy << 1097 G4double y = trueStepLength/currentRadLength;
1080 G4double y = trueStepLength / currentRadLen <<
1081 1098
1082 if(particle == positron) 1099 if(particle == positron)
1083 { 1100 {
1084 static constexpr G4double xl = 0.6; << 1101 static const G4double xl= 0.6;
1085 static constexpr G4double xh = 0.9; << 1102 static const G4double xh= 0.9;
1086 static constexpr G4double e = 113.0; << 1103 static const G4double e = 113.0;
1087 G4double corr; 1104 G4double corr;
1088 1105
1089 G4double tau = std::sqrt(currentKinEnergy << 1106 G4double tau = std::sqrt(currentKinEnergy*KineticEnergy)/mass;
1090 G4double x = std::sqrt(tau * (tau + 2.) << 1107 G4double x = std::sqrt(tau*(tau+2.)/((tau+1.)*(tau+1.)));
1091 G4double a = 0.994 - 4.08e-3 * Zeff; << 1108 G4double a = 0.994-4.08e-3*Zeff;
1092 G4double b = 7.16 + (52.6 + 365. / Zeff << 1109 G4double b = 7.16+(52.6+365./Zeff)/Zeff;
1093 G4double c = 1.000 - 4.47e-3 * Zeff; << 1110 G4double c = 1.000-4.47e-3*Zeff;
1094 G4double d = 1.21e-3 * Zeff; << 1111 G4double d = 1.21e-3*Zeff;
1095 if(x < xl) << 1112 if(x < xl) {
1096 { << 1113 corr = a*(1.-G4Exp(-b*x));
1097 corr = a * (1. - G4Exp(-b * x)); << 1114 } else if(x > xh) {
1098 } << 1115 corr = c+d*G4Exp(e*(x-1.));
1099 else if(x > xh) << 1116 } else {
1100 { << 1117 G4double yl = a*(1.-G4Exp(-b*xl));
1101 corr = c + d * G4Exp(e * (x - 1.)); << 1118 G4double yh = c+d*G4Exp(e*(xh-1.));
1102 } << 1119 G4double y0 = (yh-yl)/(xh-xl);
1103 else << 1120 G4double y1 = yl-y0*xl;
1104 { << 1121 corr = y0*x+y1;
1105 G4double yl = a * (1. - G4Exp(-b * xl)) << 1122 }
1106 G4double yh = c + d * G4Exp(e * (xh - 1 << 1123 //==================================================================
1107 G4double y0 = (yh - yl) / (xh - xl); << 1124 y *= corr*(1.+Zeff*(1.84035e-4*Zeff-1.86427e-2)+0.41125);
1108 G4double y1 = yl - y0 * xl; << 1125 }
1109 corr = y0 * x + y1; << 1126
1110 } << 1127 static const G4double c_highland = 13.6*CLHEP::MeV;
1111 y *= corr * (1. + Zeff * (1.84035e-4 * Ze << 1128 G4double theta0 = c_highland*std::abs(charge)*std::sqrt(y)*invbetacp;
1112 } << 1129
1113 <<
1114 static constexpr G4double c_highland = 13.6 <<
1115 G4double theta0 = c_highland * std::abs(cha <<
1116 <<
1117 // correction factor from e- scattering dat 1130 // correction factor from e- scattering data
1118 theta0 *= (coeffth1 + coeffth2 * G4Log(y)); << 1131 theta0 *= (coeffth1+coeffth2*G4Log(y));
1119 return theta0; 1132 return theta0;
1120 } 1133 }
1121 1134
1122 //....oooOO0OOooo........oooOO0OOooo........o 1135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 1136
1123 void G4UrbanAdjointMscModel::SampleDisplaceme 1137 void G4UrbanAdjointMscModel::SampleDisplacement(G4double sth, G4double phi)
1124 { 1138 {
1125 G4double rmax = << 1139 G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1126 sqrt((tPathLength - zPathLength) * (tPath <<
1127 <<
1128 static constexpr G4double third = 1. / 3.; <<
1129 G4double r = rmax * G4Exp(G4Log(rndmEngineM <<
1130 <<
1131 if(r > 0.) <<
1132 { <<
1133 static constexpr G4double kappa = 2.5; <<
1134 static constexpr G4double kappami1 = 1.5; <<
1135 1140
>> 1141 static const G4double third = 1./3.;
>> 1142 G4double r = rmax*G4Exp(G4Log(rndmEngineMod->flat())*third);
>> 1143 /*
>> 1144 G4cout << "G4UrbanAdjointMscModel::SampleSecondaries: e(MeV)= " << kineticEnergy
>> 1145 << " sinTheta= " << sth << " r(mm)= " << r
>> 1146 << " trueStep(mm)= " << tPathLength
>> 1147 << " geomStep(mm)= " << zPathLength
>> 1148 << G4endl;
>> 1149 */
>> 1150
>> 1151 if(r > 0.) {
>> 1152 static const G4double kappa = 2.5;
>> 1153 static const G4double kappami1 = 1.5;
>> 1154
1136 G4double latcorr = 0.; 1155 G4double latcorr = 0.;
1137 if((currentTau >= tausmall) && !insideski << 1156 if((currentTau >= tausmall) && !insideskin) {
1138 { << 1157 if(currentTau < taulim) {
1139 if(currentTau < taulim) << 1158 latcorr = lambdaeff*kappa*currentTau*currentTau*
1140 { << 1159 (1.-(kappa+1.)*currentTau*third)*third;
1141 latcorr = lambdaeff * kappa * current << 1160
1142 (1. - (kappa + 1.) * curren << 1161 } else {
1143 } << 1162 G4double etau = 0.;
1144 else << 1163 if(currentTau < taubig) { etau = G4Exp(-currentTau); }
1145 { << 1164 latcorr = -kappa*currentTau;
1146 G4double etau = 0.; << 1165 latcorr = G4Exp(latcorr)/kappami1;
1147 if(currentTau < taubig) << 1166 latcorr += 1.-kappa*etau/kappami1 ;
1148 { << 1167 latcorr *= 2.*lambdaeff*third;
1149 etau = G4Exp(-currentTau); <<
1150 } <<
1151 latcorr = -kappa * currentTau; <<
1152 latcorr = G4Exp(latcorr) / kappami1; <<
1153 latcorr += 1. - kappa * etau / kappam <<
1154 latcorr *= 2. * lambdaeff * third; <<
1155 } 1168 }
1156 } 1169 }
1157 latcorr = std::min(latcorr, r); 1170 latcorr = std::min(latcorr, r);
1158 1171
1159 // sample direction of lateral displaceme 1172 // sample direction of lateral displacement
1160 // compute it from the lateral correlatio 1173 // compute it from the lateral correlation
1161 G4double Phi = 0.; 1174 G4double Phi = 0.;
1162 if(std::abs(r * sth) < latcorr) << 1175 if(std::abs(r*sth) < latcorr) {
1163 { << 1176 Phi = twopi*rndmEngineMod->flat();
1164 Phi = twopi * rndmEngineMod->flat(); << 1177
1165 } << 1178 } else {
1166 else << 1179 //G4cout << "latcorr= " << latcorr << " r*sth= " << r*sth
1167 { << 1180 // << " ratio= " << latcorr/(r*sth) << G4endl;
1168 G4double psi = std::acos(latcorr / (r * << 1181 G4double psi = std::acos(latcorr/(r*sth));
1169 if(rndmEngineMod->flat() < 0.5) << 1182 if(rndmEngineMod->flat() < 0.5) {
1170 { << 1183 Phi = phi+psi;
1171 Phi = phi + psi; << 1184 } else {
1172 } << 1185 Phi = phi-psi;
1173 else <<
1174 { <<
1175 Phi = phi - psi; <<
1176 } 1186 }
1177 } 1187 }
1178 fDisplacement.set(r * std::cos(Phi), r * << 1188 fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1179 } 1189 }
1180 } 1190 }
1181 1191
1182 //....oooOO0OOooo........oooOO0OOooo........o 1192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1183 void G4UrbanAdjointMscModel::SampleDisplaceme << 1193
>> 1194 void G4UrbanAdjointMscModel::SampleDisplacementNew(G4double , G4double phi)
1184 { 1195 {
1185 // sample displacement r << 1196 //sample displacement r
1186 1197
1187 G4double rmax = << 1198 G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1188 sqrt((tPathLength - zPathLength) * (tPath <<
1189 // u = (r/rmax)**2 , v=1-u 1199 // u = (r/rmax)**2 , v=1-u
1190 // paramerization from ss simulation 1200 // paramerization from ss simulation
1191 // f(u) = p0*exp(p1*log(v)-p2*v)+v*(p3+p4*v << 1201 // f(u) = p0*exp(p1*log(v)-p2*v)+v*(p3+p4*v)
1192 G4double u, v, rej; << 1202 G4double u ,v , rej;
1193 G4int count = 0; 1203 G4int count = 0;
1194 1204
1195 static constexpr G4double reps = 1.e-6; << 1205 static const G4double reps = 1.e-6;
1196 static constexpr G4double rp0 = 2.2747e+4; << 1206 static const G4double rp0 = 2.2747e+4;
1197 static constexpr G4double rp1 = 4.5980e+0; << 1207 static const G4double rp1 = 4.5980e+0;
1198 static constexpr G4double rp2 = 1.5580e+1; << 1208 static const G4double rp2 = 1.5580e+1;
1199 static constexpr G4double rp3 = 7.1287e-1; << 1209 static const G4double rp3 = 7.1287e-1;
1200 static constexpr G4double rp4 = -5.7069e-1 << 1210 static const G4double rp4 =-5.7069e-1;
1201 << 1211
1202 do << 1212 do {
1203 { << 1213 u = reps+(1.-2.*reps)*rndmEngineMod->flat();
1204 u = reps + (1. - 2. * reps) * rndmEngin << 1214 v = 1.-u ;
1205 v = 1. - u; << 1215 rej = rp0*G4Exp(rp1*G4Log(v)-rp2*v) + v*(rp3+rp4*v);
1206 rej = rp0 * G4Exp(rp1 * G4Log(v) - rp2 * << 1216 }
1207 } while(rndmEngineMod->flat() > rej && ++co << 1217 // Loop checking, 15-Sept-2015, Vladimir Ivanchenko
1208 G4double r = rmax * sqrt(u); << 1218 while (rndmEngineMod->flat() > rej && ++count < 1000);
>> 1219 G4double r = rmax*sqrt(u);
1209 1220
1210 if(r > 0.) 1221 if(r > 0.)
1211 { 1222 {
1212 // sample Phi using lateral correlation 1223 // sample Phi using lateral correlation
1213 // v = Phi-phi = acos(latcorr/(r*sth)) 1224 // v = Phi-phi = acos(latcorr/(r*sth))
1214 // v has a universal distribution which c 1225 // v has a universal distribution which can be parametrized from ss
1215 // simulation as 1226 // simulation as
1216 // f(v) = 1.49e-2*exp(-v**2/(2*0.320))+2. 1227 // f(v) = 1.49e-2*exp(-v**2/(2*0.320))+2.50e-2*exp(-31.0*log(1.+6.30e-2*v))+
1217 // 1.96e-5*exp(8.42e-1*log(1.+1.45 1228 // 1.96e-5*exp(8.42e-1*log(1.+1.45e1*v))
1218 static constexpr G4double probv1 = 0.3055 << 1229 static const G4double probv1 = 0.305533;
1219 static constexpr G4double probv2 = 0.9551 << 1230 static const G4double probv2 = 0.955176;
1220 static constexpr G4double vhigh = 3.15; << 1231 static const G4double vhigh = 3.15;
1221 static const G4double w2v = 1. / G4Exp(30 << 1232 static const G4double w2v = 1./G4Exp(30.*G4Log(1. + 6.30e-2*vhigh));
1222 static const G4double w3v = 1. / G4Exp(-1 << 1233 static const G4double w3v = 1./G4Exp(-1.842*G4Log(1. + 1.45e1*vhigh));
1223 1234
1224 G4double Phi; 1235 G4double Phi;
1225 G4double random = rndmEngineMod->flat(); 1236 G4double random = rndmEngineMod->flat();
1226 if(random < probv1) << 1237 if(random < probv1) {
1227 { << 1238 do {
1228 do << 1239 v = G4RandGauss::shoot(rndmEngineMod,0.,0.320);
1229 { <<
1230 v = G4RandGauss::shoot(rndmEngineMod, <<
1231 } while(std::abs(v) >= vhigh); <<
1232 Phi = phi + v; <<
1233 } <<
1234 else <<
1235 { <<
1236 if(random < probv2) <<
1237 { <<
1238 v = (-1. + <<
1239 1. / G4Exp(G4Log(1. - rndmEngine <<
1240 6.30e-2; <<
1241 } 1240 }
1242 else << 1241 // Loop checking, 15-Sept-2015, Vladimir Ivanchenko
1243 { << 1242 while (std::abs(v) >= vhigh);
1244 v = (-1. + 1. / G4Exp(G4Log(1. - rndm << 1243 Phi = phi + v;
1245 -1.842)) / << 1244
1246 1.45e1; << 1245 } else {
>> 1246
>> 1247 if(random < probv2) {
>> 1248 v = (-1.+1./G4Exp(G4Log(1.-rndmEngineMod->flat()*(1.-w2v))/30.))/6.30e-2;
>> 1249 } else {
>> 1250 v = (-1.+1./G4Exp(G4Log(1.-rndmEngineMod->flat()*(1.-w3v))/-1.842))/1.45e1;
1247 } 1251 }
1248 1252
1249 random = rndmEngineMod->flat(); 1253 random = rndmEngineMod->flat();
1250 if(random < 0.5) << 1254 if(random < 0.5) { Phi = phi+v; }
1251 { << 1255 else { Phi = phi-v; }
1252 Phi = phi + v; <<
1253 } <<
1254 else <<
1255 { <<
1256 Phi = phi - v; <<
1257 } <<
1258 } 1256 }
1259 fDisplacement.set(r * std::cos(Phi), r * << 1257 fDisplacement.set(r*std::cos(Phi),r*std::sin(Phi),0.0);
1260 } 1258 }
1261 } 1259 }
1262 1260
1263 //....oooOO0OOooo........oooOO0OOooo........o 1261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1264 1262