| 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: G4WentzelVIModel.cc 88979 2015-03-17 10:10:21Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class file 30 // GEANT4 Class file
30 // 31 //
31 // 32 //
32 // File name: G4WentzelVIModel 33 // File name: G4WentzelVIModel
33 // 34 //
34 // Author: V.Ivanchenko 35 // Author: V.Ivanchenko
35 // 36 //
36 // Creation date: 09.04.2008 from G4MuMscModel 37 // Creation date: 09.04.2008 from G4MuMscModel
37 // 38 //
38 // Modifications: 39 // Modifications:
39 // 27-05-2010 V.Ivanchenko added G4WentzelOKan 40 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
40 // compute cross sections and sam 41 // compute cross sections and sample scattering angle
41 // 42 //
42 // 43 //
43 // Class Description: 44 // Class Description:
44 // 45 //
45 // Implementation of the model of multiple sca 46 // Implementation of the model of multiple scattering based on
46 // G.Wentzel, Z. Phys. 40 (1927) 590. 47 // G.Wentzel, Z. Phys. 40 (1927) 590.
47 // H.W.Lewis, Phys Rev 78 (1950) 526. 48 // H.W.Lewis, Phys Rev 78 (1950) 526.
48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 49 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
49 // L.Urban, CERN-OPEN-2006-077. 50 // L.Urban, CERN-OPEN-2006-077.
50 51
51 // ------------------------------------------- 52 // -------------------------------------------------------------------
52 // 53 //
53 54
54 //....oooOO0OOooo........oooOO0OOooo........oo 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55 //....oooOO0OOooo........oooOO0OOooo........oo 56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 57
57 #include "G4WentzelVIModel.hh" 58 #include "G4WentzelVIModel.hh"
58 #include "G4PhysicalConstants.hh" 59 #include "G4PhysicalConstants.hh"
59 #include "G4SystemOfUnits.hh" 60 #include "G4SystemOfUnits.hh"
60 #include "Randomize.hh" 61 #include "Randomize.hh"
61 #include "G4ParticleChangeForMSC.hh" 62 #include "G4ParticleChangeForMSC.hh"
62 #include "G4PhysicsTableHelper.hh" 63 #include "G4PhysicsTableHelper.hh"
63 #include "G4ElementVector.hh" 64 #include "G4ElementVector.hh"
64 #include "G4ProductionCutsTable.hh" 65 #include "G4ProductionCutsTable.hh"
65 #include "G4EmParameters.hh" 66 #include "G4EmParameters.hh"
66 #include "G4Log.hh" 67 #include "G4Log.hh"
67 #include "G4Exp.hh" 68 #include "G4Exp.hh"
68 69
69 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
70 71
71 using namespace std; 72 using namespace std;
72 73
73 const G4double invsqrt12 = 1./std::sqrt(12.); <<
74 const G4double numlimit = 0.1; <<
75 const G4int minNCollisions = 10; 74 const G4int minNCollisions = 10;
76 75
77 G4WentzelVIModel::G4WentzelVIModel(G4bool comb << 76 G4WentzelVIModel::G4WentzelVIModel(G4bool combined, const G4String& nam) :
78 : G4VMscModel(nam), << 77 G4VMscModel(nam),
79 singleScatteringMode(false), << 78 ssFactor(1.05),
80 isCombined(comb), << 79 invssFactor(1.0),
81 useSecondMoment(false) << 80 currentCouple(0),
>> 81 inside(false),
>> 82 singleScatteringMode(false),
>> 83 cosThetaMin(1.0),
>> 84 cosThetaMax(-1.0),
>> 85 fSecondMoments(0),
>> 86 idx2(0),
>> 87 numlimit(0.1),
>> 88 isCombined(combined),
>> 89 useSecondMoment(false)
82 { 90 {
83 tlimitminfix = 1.e-6*CLHEP::mm; <<
84 lowEnergyLimit = 1.0*CLHEP::eV; <<
85 SetSingleScatteringFactor(1.25); 91 SetSingleScatteringFactor(1.25);
86 wokvi = new G4WentzelOKandVIxSection(isCombi << 92 invsqrt12 = 1./sqrt(12.);
>> 93 tlimitminfix = 1.e-6*mm;
>> 94 lowEnergyLimit = 1.0*eV;
>> 95 particle = 0;
>> 96 nelments = 5;
>> 97 xsecn.resize(nelments);
>> 98 prob.resize(nelments);
>> 99 wokvi = new G4WentzelOKandVIxSection(combined);
>> 100 fixedCut = -1.0;
>> 101
>> 102 preKinEnergy = effKinEnergy = tPathLength = zPathLength = lambdaeff
>> 103 = currentRange = xtsec = cosTetMaxNuc = 0.0;
>> 104 currentMaterialIndex = 0;
>> 105
>> 106 fParticleChange = 0;
>> 107 currentCuts = 0;
>> 108 currentMaterial = 0;
87 } 109 }
88 110
89 //....oooOO0OOooo........oooOO0OOooo........oo 111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
90 112
91 G4WentzelVIModel::~G4WentzelVIModel() 113 G4WentzelVIModel::~G4WentzelVIModel()
92 { 114 {
93 delete wokvi; 115 delete wokvi;
94 if(IsMaster()) { << 116 if(fSecondMoments && IsMaster()) {
95 delete fSecondMoments; 117 delete fSecondMoments;
96 fSecondMoments = nullptr; << 118 fSecondMoments = 0;
97 } 119 }
98 } 120 }
99 121
100 //....oooOO0OOooo........oooOO0OOooo........oo 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 123
102 void G4WentzelVIModel::Initialise(const G4Part 124 void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
103 const G4Data << 125 const G4DataVector& cuts)
104 { 126 {
105 // reset parameters 127 // reset parameters
106 SetupParticle(p); 128 SetupParticle(p);
107 InitialiseParameters(p); <<
108 currentRange = 0.0; 129 currentRange = 0.0;
109 130
110 if(isCombined) { 131 if(isCombined) {
111 G4double tet = PolarAngleLimit(); 132 G4double tet = PolarAngleLimit();
112 if(tet <= 0.0) { cosThetaMax = 1 133 if(tet <= 0.0) { cosThetaMax = 1.0; }
113 else if(tet < CLHEP::pi) { cosThetaMax = c 134 else if(tet < CLHEP::pi) { cosThetaMax = cos(tet); }
114 } 135 }
115 //G4cout << "G4WentzelVIModel::Initialise " <<
116 // << " " << this << " " << wokvi << G4end <<
117 136
>> 137 //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
118 wokvi->Initialise(p, cosThetaMax); 138 wokvi->Initialise(p, cosThetaMax);
119 /* << 139 /*
120 G4cout << "G4WentzelVIModel: " << particle-> 140 G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
121 << " 1-cos(ThetaLimit)= " << 1 - cos 141 << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
122 << " SingScatFactor= " << ssFactor << 142 << " SingScatFactor= " << ssFactor
123 << G4endl; << 143 << G4endl;
124 */ 144 */
125 currentCuts = &cuts; 145 currentCuts = &cuts;
126 146
127 // set values of some data members 147 // set values of some data members
128 fParticleChange = GetParticleChangeForMSC(p) 148 fParticleChange = GetParticleChangeForMSC(p);
129 149
130 // Access to materials <<
131 const G4ProductionCutsTable* theCoupleTable <<
132 G4ProductionCutsTable::GetProductionCutsTa <<
133 G4int numOfCouples = (G4int)theCoupleTable-> <<
134 nelments = 0; <<
135 for(G4int i=0; i<numOfCouples; ++i) { <<
136 G4int nelm = (G4int)theCoupleTable->GetMat <<
137 nelments = std::max(nelments, nelm); <<
138 } <<
139 xsecn.resize(nelments); <<
140 prob.resize(nelments); <<
141 <<
142 // build second moment table only if transpo 150 // build second moment table only if transport table is build
143 G4PhysicsTable* table = GetCrossSectionTable 151 G4PhysicsTable* table = GetCrossSectionTable();
144 if(useSecondMoment && IsMaster() && nullptr << 152 if(useSecondMoment && IsMaster() && table) {
145 153
146 //G4cout << "### G4WentzelVIModel::Initial 154 //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table "
147 // << table << G4endl; << 155 // << table << G4endl;
148 fSecondMoments = 156 fSecondMoments =
149 G4PhysicsTableHelper::PreparePhysicsTabl 157 G4PhysicsTableHelper::PreparePhysicsTable(fSecondMoments);
>> 158 // Access to materials
>> 159 const G4ProductionCutsTable* theCoupleTable =
>> 160 G4ProductionCutsTable::GetProductionCutsTable();
>> 161 size_t numOfCouples = theCoupleTable->GetTableSize();
150 162
151 G4bool splineFlag = true; 163 G4bool splineFlag = true;
152 G4PhysicsVector* aVector = nullptr; << 164 G4PhysicsVector* aVector = 0;
153 G4PhysicsVector* bVector = nullptr; << 165 G4PhysicsVector* bVector = 0;
154 G4double emin = std::max(LowEnergyLimit(), 166 G4double emin = std::max(LowEnergyLimit(), LowEnergyActivationLimit());
155 G4double emax = std::min(HighEnergyLimit() 167 G4double emax = std::min(HighEnergyLimit(), HighEnergyActivationLimit());
156 if(emin < emax) { 168 if(emin < emax) {
157 std::size_t n = G4EmParameters::Instance << 169 size_t n = G4EmParameters::Instance()->NumberOfBinsPerDecade()
158 *G4lrint(std::log10(emax/emin)); << 170 *G4lrint(std::log10(emax/emin));
159 if(n < 3) { n = 3; } 171 if(n < 3) { n = 3; }
160 172
161 for(G4int i=0; i<numOfCouples; ++i) { << 173 for(size_t i=0; i<numOfCouples; ++i) {
162 174
163 //G4cout<< "i= " << i << " Flag= " << << 175 //G4cout<< "i= " << i << " Flag= " << fSecondMoments->GetFlag(i)
164 // << G4endl; << 176 // << G4endl;
165 if(fSecondMoments->GetFlag(i)) { << 177 if(fSecondMoments->GetFlag(i)) {
166 DefineMaterial(theCoupleTable->GetMa << 178 DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
167 179
168 delete (*fSecondMoments)[i]; << 180 delete (*fSecondMoments)[i];
169 if(nullptr == aVector) { << 181 if(!aVector) {
170 aVector = new G4PhysicsLogVector(e << 182 aVector = new G4PhysicsLogVector(emin, emax, n);
171 bVector = aVector; << 183 bVector = aVector;
172 } else { << 184 } else {
173 bVector = new G4PhysicsVector(*aVe << 185 bVector = new G4PhysicsVector(*aVector);
174 } << 186 }
175 for(std::size_t j=0; j<n; ++j) { << 187 for(size_t j=0; j<n; ++j) {
176 G4double e = bVector->Energy(j); << 188 G4double e = bVector->Energy(j);
177 bVector->PutValue(j, ComputeSecond << 189 bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e);
178 } << 190 }
179 if(splineFlag) { bVector->FillSecond << 191 if(splineFlag) { bVector->FillSecondDerivatives(); }
180 (*fSecondMoments)[i] = bVector; << 192 (*fSecondMoments)[i] = bVector;
181 } << 193 }
182 } 194 }
183 } 195 }
184 //G4cout << *fSecondMoments << G4endl; 196 //G4cout << *fSecondMoments << G4endl;
185 } 197 }
186 } 198 }
187 199
188 //....oooOO0OOooo........oooOO0OOooo........oo 200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 201
190 void G4WentzelVIModel::InitialiseLocal(const G 202 void G4WentzelVIModel::InitialiseLocal(const G4ParticleDefinition*,
191 G4VEmMo << 203 G4VEmModel* masterModel)
192 { 204 {
193 fSecondMoments = static_cast<G4WentzelVIMode 205 fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel)
194 ->GetSecondMomentTable(); 206 ->GetSecondMomentTable();
195 } 207 }
196 208
197 //....oooOO0OOooo........oooOO0OOooo........oo <<
198 <<
199 void G4WentzelVIModel::DefineMaterial(const G4 <<
200 { <<
201 if(cup != currentCouple) { <<
202 currentCouple = cup; <<
203 SetCurrentCouple(cup); <<
204 currentMaterial = cup->GetMaterial(); <<
205 currentMaterialIndex = currentCouple->GetI <<
206 } <<
207 } <<
208 <<
209 //....oooOO0OOooo........oooOO0OOooo........oo 209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 210
211 G4double G4WentzelVIModel::ComputeCrossSection 211 G4double G4WentzelVIModel::ComputeCrossSectionPerAtom(
212 const G4ParticleD 212 const G4ParticleDefinition* p,
213 G4double kinEnerg << 213 G4double kinEnergy,
214 G4double Z, G4dou << 214 G4double Z, G4double,
215 G4double cutEnerg << 215 G4double cutEnergy, G4double)
216 { 216 {
217 G4double cross = 0.0; 217 G4double cross = 0.0;
218 SetupParticle(p); << 218 if(p != particle) { SetupParticle(p); }
219 if(kinEnergy < lowEnergyLimit) { return cros 219 if(kinEnergy < lowEnergyLimit) { return cross; }
220 if(nullptr == CurrentCouple()) { << 220 if(!CurrentCouple()) {
221 G4Exception("G4WentzelVIModel::ComputeCros 221 G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
222 FatalException, " G4MaterialCu << 222 FatalException, " G4MaterialCutsCouple is not defined");
223 return 0.0; 223 return 0.0;
224 } 224 }
225 DefineMaterial(CurrentCouple()); 225 DefineMaterial(CurrentCouple());
226 cosTetMaxNuc = wokvi->SetupKinematic(kinEner 226 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
227 if(cosTetMaxNuc < 1.0) { 227 if(cosTetMaxNuc < 1.0) {
228 G4double cut = (0.0 < fixedCut) ? fixedCu << 228 G4double cut = cutEnergy;
>> 229 if(fixedCut > 0.0) { cut = fixedCut; }
229 G4double cost = wokvi->SetupTarget(G4lrint 230 G4double cost = wokvi->SetupTarget(G4lrint(Z), cut);
230 cross = wokvi->ComputeTransportCrossSectio 231 cross = wokvi->ComputeTransportCrossSectionPerAtom(cost);
231 /* 232 /*
232 if(p->GetParticleName() == "e-") 233 if(p->GetParticleName() == "e-")
233 G4cout << "G4WentzelVIModel::CS: Z= " << G 234 G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy
234 << " 1-cosN= " << 1 - cosTetMaxNuc << 235 << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
235 << " " << particle->GetParticleName << 236 << " " << particle->GetParticleName() << G4endl;
236 */ 237 */
237 } 238 }
238 return cross; 239 return cross;
239 } 240 }
240 241
241 //....oooOO0OOooo........oooOO0OOooo........oo 242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 243
243 void G4WentzelVIModel::StartTracking(G4Track* 244 void G4WentzelVIModel::StartTracking(G4Track* track)
244 { 245 {
245 /* << 246 SetupParticle(track->GetDynamicParticle()->GetDefinition());
246 G4cout << "G4WentzelVIModel::StartTracking " << 247 inside = false;
247 << track->GetParticleDefinition()->GetParti << 248 G4VEmModel::StartTracking(track);
248 << " workvi: " << wokvi << G4endl; <<
249 */ <<
250 SetupParticle(track->GetParticleDefinition() <<
251 } 249 }
252 250
253 //....oooOO0OOooo........oooOO0OOooo........oo 251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 252
255 G4double G4WentzelVIModel::ComputeTruePathLeng 253 G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
256 const G4Track& tr 254 const G4Track& track,
257 G4double& current << 255 G4double& currentMinimalStep)
258 { 256 {
259 G4double tlimit = currentMinimalStep; 257 G4double tlimit = currentMinimalStep;
260 const G4DynamicParticle* dp = track.GetDynam 258 const G4DynamicParticle* dp = track.GetDynamicParticle();
261 const G4StepPoint* sp = track.GetStep()->Get << 259 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
262 G4StepStatus stepStatus = sp->GetStepStatus( 260 G4StepStatus stepStatus = sp->GetStepStatus();
263 singleScatteringMode = false; 261 singleScatteringMode = false;
264 262
265 //G4cout << "G4WentzelVIModel::ComputeTruePa 263 //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
266 // << stepStatus << " " << track.Ge << 264 // << stepStatus << " " << track.GetDefinition()->GetParticleName()
267 // << G4endl; << 265 // << G4endl;
268 266
269 // initialisation for each step, lambda may 267 // initialisation for each step, lambda may be computed from scratch
270 preKinEnergy = dp->GetKineticEnergy(); 268 preKinEnergy = dp->GetKineticEnergy();
271 effKinEnergy = preKinEnergy; 269 effKinEnergy = preKinEnergy;
272 DefineMaterial(track.GetMaterialCutsCouple() 270 DefineMaterial(track.GetMaterialCutsCouple());
273 const G4double logPreKinEnergy = dp->GetLogK << 271 lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
274 lambdaeff = GetTransportMeanFreePath(particl << 272 currentRange = GetRange(particle,preKinEnergy,currentCouple);
275 currentRange = GetRange(particle,preKinEnerg <<
276 cosTetMaxNuc = wokvi->SetupKinematic(preKinE 273 cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
277 274
278 //G4cout << "lambdaeff= " << lambdaeff << " 275 //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
279 // << " tlimit= " << tlimit << " 1-cost= " < << 276 // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
280 277
281 // extra check for abnormal situation 278 // extra check for abnormal situation
282 // this check needed to run MSC with eIoni a 279 // this check needed to run MSC with eIoni and eBrem inactivated
283 if(tlimit > currentRange) { tlimit = current 280 if(tlimit > currentRange) { tlimit = currentRange; }
284 281
285 // stop here if small range particle 282 // stop here if small range particle
286 if(tlimit < tlimitminfix) { << 283 if(inside || tlimit < tlimitminfix) {
287 return ConvertTrueToGeom(tlimit, currentMi 284 return ConvertTrueToGeom(tlimit, currentMinimalStep);
288 } 285 }
289 286
290 // pre step 287 // pre step
291 G4double presafety = sp->GetSafety(); 288 G4double presafety = sp->GetSafety();
292 // far from geometry boundary 289 // far from geometry boundary
293 if(currentRange < presafety) { 290 if(currentRange < presafety) {
>> 291 inside = true;
294 return ConvertTrueToGeom(tlimit, currentMi 292 return ConvertTrueToGeom(tlimit, currentMinimalStep);
295 } 293 }
296 294
297 // compute presafety again if presafety <= 0 295 // compute presafety again if presafety <= 0 and no boundary
298 // i.e. when it is needed for optimization p 296 // i.e. when it is needed for optimization purposes
299 if(stepStatus != fGeomBoundary && presafety 297 if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
300 presafety = ComputeSafety(sp->GetPosition( 298 presafety = ComputeSafety(sp->GetPosition(), tlimit);
301 if(currentRange < presafety) { 299 if(currentRange < presafety) {
>> 300 inside = true;
302 return ConvertTrueToGeom(tlimit, current 301 return ConvertTrueToGeom(tlimit, currentMinimalStep);
303 } 302 }
304 } 303 }
305 /* 304 /*
306 G4cout << "e(MeV)= " << preKinEnergy/MeV 305 G4cout << "e(MeV)= " << preKinEnergy/MeV
307 << " " << particle->GetParticleName( << 306 << " " << particle->GetParticleName()
308 << " CurLimit(mm)= " << tlimit/mm <<" << 307 << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
309 << " R(mm)= " <<currentRange/mm << 308 << " R(mm)= " <<currentRange/mm
310 << " L0(mm^-1)= " << lambdaeff*mm << 309 << " L0(mm^-1)= " << lambdaeff*mm
311 << G4endl; << 310 << G4endl;
312 */ 311 */
313 // natural limit for high energy 312 // natural limit for high energy
314 G4double rlimit = std::max(facrange*currentR 313 G4double rlimit = std::max(facrange*currentRange,
315 (1.0 - cosTetMaxN << 314 (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor);
316 315
317 // low-energy e- 316 // low-energy e-
318 if(cosThetaMax > cosTetMaxNuc) { 317 if(cosThetaMax > cosTetMaxNuc) {
319 rlimit = std::min(rlimit, facsafety*presaf 318 rlimit = std::min(rlimit, facsafety*presafety);
320 } 319 }
321 320
322 // cut correction 321 // cut correction
323 G4double rcut = currentCouple->GetProduction 322 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
324 //G4cout << "rcut= " << rcut << " rlimit= " 323 //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
325 // << presafety << " 1-cosThetaMax= " <<1-co 324 // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
326 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc < 325 //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
327 if(rcut > rlimit) { rlimit = std::min(rlimit << 326 if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
328 327
329 tlimit = std::min(tlimit, rlimit); 328 tlimit = std::min(tlimit, rlimit);
330 tlimit = std::max(tlimit, tlimitminfix); 329 tlimit = std::max(tlimit, tlimitminfix);
331 330
332 // step limit in infinite media 331 // step limit in infinite media
333 tlimit = std::min(tlimit, 50*currentMaterial 332 tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
334 333
335 //compute geomlimit and force few steps with 334 //compute geomlimit and force few steps within a volume
336 if (steppingAlgorithm == fUseDistanceToBound 335 if (steppingAlgorithm == fUseDistanceToBoundary
337 && stepStatus == fGeomBoundary) { 336 && stepStatus == fGeomBoundary) {
338 337
339 G4double geomlimit = ComputeGeomLimit(trac 338 G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
340 tlimit = std::min(tlimit, geomlimit/facgeo 339 tlimit = std::min(tlimit, geomlimit/facgeom);
341 } 340 }
342 /* 341 /*
343 G4cout << particle->GetParticleName() << " e 342 G4cout << particle->GetParticleName() << " e= " << preKinEnergy
344 << " L0= " << lambdaeff << " R= " << << 343 << " L0= " << lambdaeff << " R= " << currentRange
345 << " tlimit= " << tlimit << 344 << " tlimit= " << tlimit
346 << " currentMinimalStep= " << curre << 345 << " currentMinimalStep= " << currentMinimalStep << G4endl;
347 */ 346 */
348 return ConvertTrueToGeom(tlimit, currentMini 347 return ConvertTrueToGeom(tlimit, currentMinimalStep);
349 } 348 }
350 349
351 //....oooOO0OOooo........oooOO0OOooo........oo 350 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352 351
353 G4double G4WentzelVIModel::ComputeGeomPathLeng 352 G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
354 { 353 {
355 zPathLength = tPathLength = truelength; 354 zPathLength = tPathLength = truelength;
356 355
357 // small step use only single scattering 356 // small step use only single scattering
358 cosThetaMin = 1.0; 357 cosThetaMin = 1.0;
359 ComputeTransportXSectionPerVolume(cosThetaMi 358 ComputeTransportXSectionPerVolume(cosThetaMin);
360 //G4cout << "xtsec= " << xtsec << " Nav= " 359 //G4cout << "xtsec= " << xtsec << " Nav= "
361 // << zPathLength*xtsec << G4endl; << 360 // << zPathLength*xtsec << G4endl;
362 if(0.0 >= lambdaeff || G4int(zPathLength*xts 361 if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
363 singleScatteringMode = true; 362 singleScatteringMode = true;
364 lambdaeff = DBL_MAX; 363 lambdaeff = DBL_MAX;
365 364
366 } else { 365 } else {
367 //G4cout << "ComputeGeomPathLength: tLengt 366 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
368 // << " Leff= " << lambdaeff << << 367 // << " Leff= " << lambdaeff << G4endl;
369 // small step 368 // small step
370 if(tPathLength < numlimit*lambdaeff) { 369 if(tPathLength < numlimit*lambdaeff) {
371 G4double tau = tPathLength/lambdaeff; 370 G4double tau = tPathLength/lambdaeff;
372 zPathLength *= (1.0 - 0.5*tau + tau*tau/ 371 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
373 372
374 // medium step 373 // medium step
375 } else { 374 } else {
376 G4double e1 = 0.0; 375 G4double e1 = 0.0;
377 if(currentRange > tPathLength) { 376 if(currentRange > tPathLength) {
378 e1 = GetEnergy(particle,currentRange-t << 377 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
379 } 378 }
380 effKinEnergy = 0.5*(e1 + preKinEnergy); 379 effKinEnergy = 0.5*(e1 + preKinEnergy);
381 cosTetMaxNuc = wokvi->SetupKinematic(eff 380 cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial);
382 lambdaeff = GetTransportMeanFreePath(par 381 lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy);
383 //G4cout << " tLength= "<< tPathLength<< 382 //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
384 zPathLength = lambdaeff; 383 zPathLength = lambdaeff;
385 if(tPathLength*numlimit < lambdaeff) { 384 if(tPathLength*numlimit < lambdaeff) {
386 zPathLength *= (1.0 - G4Exp(-tPathLeng << 385 zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
387 } 386 }
388 } 387 }
389 } 388 }
390 //G4cout << "Comp.geom: zLength= "<<zPathLen 389 //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
391 // << tPathLength<< " Leff= " << lam << 390 // << tPathLength<< " Leff= " << lambdaeff << G4endl;
392 return zPathLength; 391 return zPathLength;
393 } 392 }
394 393
395 //....oooOO0OOooo........oooOO0OOooo........oo 394 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396 395
397 G4double G4WentzelVIModel::ComputeTrueStepLeng 396 G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
398 { 397 {
399 // initialisation of single scattering x-sec 398 // initialisation of single scattering x-section
400 /* 399 /*
401 G4cout << "ComputeTrueStepLength: Step= " << 400 G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
402 << " geomL= " << zPathLength << 401 << " geomL= " << zPathLength
403 << " Lambda= " << lambdaeff << 402 << " Lambda= " << lambdaeff
404 << " 1-cosThetaMaxNuc= " << 1 - cos << 403 << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
405 */ 404 */
406 if(singleScatteringMode) { 405 if(singleScatteringMode) {
407 zPathLength = tPathLength = geomStepLength 406 zPathLength = tPathLength = geomStepLength;
408 407
409 } else { 408 } else {
410 409
411 // step defined by transportation 410 // step defined by transportation
412 // change both geom and true step lengths 411 // change both geom and true step lengths
413 if(geomStepLength < zPathLength) { 412 if(geomStepLength < zPathLength) {
414 413
415 // single scattering 414 // single scattering
416 if(G4int(geomStepLength*xtsec) < minNCol 415 if(G4int(geomStepLength*xtsec) < minNCollisions) {
417 zPathLength = tPathLength = geomStepLe << 416 zPathLength = tPathLength = geomStepLength;
418 lambdaeff = DBL_MAX; << 417 lambdaeff = DBL_MAX;
419 singleScatteringMode = true; << 418 singleScatteringMode = true;
420 419
421 // multiple scattering << 420 // multiple scattering
422 } else { 421 } else {
423 // small step << 422 // small step
424 if(geomStepLength < numlimit*lambdaeff << 423 if(geomStepLength < numlimit*lambdaeff) {
425 G4double tau = geomStepLength/lambda << 424 G4double tau = geomStepLength/lambdaeff;
426 tPathLength = geomStepLength*(1.0 + << 425 tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0);
427 << 426
428 // energy correction for a big step << 427 // energy correction for a big step
429 } else { << 428 } else {
430 tPathLength *= geomStepLength/zPathL << 429 tPathLength *= geomStepLength/zPathLength;
431 G4double e1 = 0.0; << 430 G4double e1 = 0.0;
432 if(currentRange > tPathLength) { << 431 if(currentRange > tPathLength) {
433 e1 = GetEnergy(particle,currentRan << 432 e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
434 } << 433 }
435 effKinEnergy = 0.5*(e1 + preKinEnerg << 434 effKinEnergy = 0.5*(e1 + preKinEnergy);
436 cosTetMaxNuc = wokvi->SetupKinematic << 435 cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial);
437 lambdaeff = GetTransportMeanFreePath << 436 lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy);
438 G4double tau = geomStepLength/lambda << 437 G4double tau = geomStepLength/lambdaeff;
439 << 438
440 if(tau < 0.999999) { tPathLength = - << 439 if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
441 else { tPathLength = c << 440 else { tPathLength = currentRange; }
442 } << 441 }
443 zPathLength = geomStepLength; << 442 zPathLength = geomStepLength;
444 } 443 }
445 } 444 }
446 } 445 }
447 // check of step length 446 // check of step length
448 // define threshold angle between single and 447 // define threshold angle between single and multiple scattering
449 if(!singleScatteringMode) { 448 if(!singleScatteringMode) {
450 cosThetaMin -= ssFactor*tPathLength/lambda 449 cosThetaMin -= ssFactor*tPathLength/lambdaeff;
451 xtsec = 0.0; 450 xtsec = 0.0;
452 451
453 // recompute transport cross section - do 452 // recompute transport cross section - do not change energy
454 // anymore - cannot be applied for big ste 453 // anymore - cannot be applied for big steps
455 if(cosThetaMin > cosTetMaxNuc) { 454 if(cosThetaMin > cosTetMaxNuc) {
456 // new computation 455 // new computation
457 G4double cross = ComputeTransportXSectio 456 G4double cross = ComputeTransportXSectionPerVolume(cosThetaMin);
458 //G4cout << "%%%% cross= " << cross << " 457 //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec
459 // << " 1-cosTMin= " << 1.0 - << 458 // << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
460 if(cross <= 0.0) { 459 if(cross <= 0.0) {
461 singleScatteringMode = true; << 460 singleScatteringMode = true;
462 tPathLength = zPathLength; << 461 tPathLength = zPathLength;
463 lambdaeff = DBL_MAX; << 462 lambdaeff = DBL_MAX;
464 cosThetaMin = 1.0; << 463 cosThetaMin = 1.0;
465 } else if(xtsec > 0.0) { 464 } else if(xtsec > 0.0) {
466 << 465
467 lambdaeff = 1./cross; << 466 lambdaeff = 1./cross;
468 G4double tau = zPathLength*cross; << 467 G4double tau = zPathLength*cross;
469 if(tau < numlimit) { << 468 if(tau < numlimit) {
470 tPathLength = zPathLength*(1.0 + 0.5 << 469 tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
471 } else if(tau < 0.999999) { << 470 } else if(tau < 0.999999) {
472 tPathLength = -lambdaeff*G4Log(1.0 - << 471 tPathLength = -lambdaeff*G4Log(1.0 - tau);
473 } else { << 472 } else {
474 tPathLength = currentRange; << 473 tPathLength = currentRange;
475 } << 474 }
476 } 475 }
477 } 476 }
478 } 477 }
479 tPathLength = std::min(tPathLength, currentR 478 tPathLength = std::min(tPathLength, currentRange);
480 /* 479 /*
481 G4cout <<"Comp.true: zLength= "<<zPathLength 480 G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
482 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0 << 481 <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
483 G4cout << particle->GetParticleName() << " 1 482 G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
484 << " 1-cosTetMaxNuc= " << 1-cosTetMax << 483 << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
485 << " e(MeV)= " << preKinEnergy/MeV << << 484 << " e(MeV)= " << preKinEnergy/MeV << " "
486 << " SSmode= " << singleScatteringMod << 485 << " SSmode= " << singleScatteringMode << G4endl;
487 */ 486 */
488 return tPathLength; 487 return tPathLength;
489 } 488 }
490 489
491 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
492 491
493 G4ThreeVector& 492 G4ThreeVector&
494 G4WentzelVIModel::SampleScattering(const G4Thr 493 G4WentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection,
495 G4double /* << 494 G4double /*safety*/)
496 { 495 {
497 fDisplacement.set(0.0,0.0,0.0); 496 fDisplacement.set(0.0,0.0,0.0);
498 //G4cout << "!##! G4WentzelVIModel::SampleSc 497 //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
499 // << particle->GetParticleName() << << 498 // << particle->GetParticleName() << G4endl;
500 499
501 // ignore scattering for zero step length an 500 // ignore scattering for zero step length and energy below the limit
502 if(preKinEnergy < lowEnergyLimit || tPathLen 501 if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
503 { return fDisplacement; } 502 { return fDisplacement; }
504 503
505 G4double invlambda = 0.0; 504 G4double invlambda = 0.0;
506 if(lambdaeff < DBL_MAX) { invlambda = 0.5/la 505 if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
507 506
508 // use average kinetic energy over the step 507 // use average kinetic energy over the step
509 G4double cut = (*currentCuts)[currentMateria 508 G4double cut = (*currentCuts)[currentMaterialIndex];
510 if(fixedCut > 0.0) { cut = fixedCut; } 509 if(fixedCut > 0.0) { cut = fixedCut; }
511 /* 510 /*
512 G4cout <<"SampleScat: E0(MeV)= "<< preKinEne 511 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
513 << " Leff= " << lambdaeff <<" sig0( << 512 << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
514 << " xmsc= " << tPathLength*invlamb << 513 << " xmsc= " << tPathLength*invlambda
515 << " safety= " << safety << G4endl; << 514 << " safety= " << safety << G4endl;
516 */ 515 */
517 // step limit due msc 516 // step limit due msc
518 G4int nMscSteps = 1; 517 G4int nMscSteps = 1;
519 G4double x0 = tPathLength; 518 G4double x0 = tPathLength;
520 G4double z0 = x0*invlambda; 519 G4double z0 = x0*invlambda;
>> 520 //G4double zzz = 0.0;
521 G4double prob2 = 0.0; 521 G4double prob2 = 0.0;
522 522
523 CLHEP::HepRandomEngine* rndmEngine = G4Rando <<
524 <<
525 // large scattering angle case - two step ap 523 // large scattering angle case - two step approach
>> 524
526 if(!singleScatteringMode) { 525 if(!singleScatteringMode) {
>> 526 static const G4double zzmin = 0.05;
527 if(useSecondMoment) { 527 if(useSecondMoment) {
528 G4double z1 = invlambda*invlambda; 528 G4double z1 = invlambda*invlambda;
529 G4double z2 = SecondMoment(particle, cur 529 G4double z2 = SecondMoment(particle, currentCouple, effKinEnergy);
530 prob2 = (z2 - z1)/(1.5*z1 - z2); 530 prob2 = (z2 - z1)/(1.5*z1 - z2);
531 } 531 }
532 // numerical limitation on step length for << 532 // if(z0 > zzmin && safety > tlimitminfix) {
533 if (z0 > 0.05 && z0 < 30.) { << 533 if(z0 > zzmin) {
534 x0 *= 0.5; 534 x0 *= 0.5;
535 z0 *= 0.5; 535 z0 *= 0.5;
536 nMscSteps = 2; 536 nMscSteps = 2;
537 G4double zzz = G4Exp(-1.0/z0); << 537 }
>> 538 //if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
>> 539 G4double zzz = 0.0;
>> 540 if(z0 > zzmin) {
>> 541 zzz = G4Exp(-1.0/z0);
538 z0 += zzz; 542 z0 += zzz;
539 prob2 *= (1.0 + zzz); << 543 prob2 *= (1 + zzz);
540 } 544 }
541 prob2 /= (1.0 + prob2); << 545 prob2 /= (1 + prob2);
542 } 546 }
543 547
544 // step limit due to single scattering 548 // step limit due to single scattering
545 G4double x1 = 2*tPathLength; 549 G4double x1 = 2*tPathLength;
546 if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->fl << 550 if(0.0 < xtsec) { x1 = -G4Log(G4UniformRand())/xtsec; }
547 551
548 // no scattering case 552 // no scattering case
549 if(singleScatteringMode && x1 > tPathLength) 553 if(singleScatteringMode && x1 > tPathLength)
550 { return fDisplacement; } 554 { return fDisplacement; }
551 555
552 const G4ElementVector* theElementVector = 556 const G4ElementVector* theElementVector =
553 currentMaterial->GetElementVector(); 557 currentMaterial->GetElementVector();
554 std::size_t nelm = currentMaterial->GetNumbe << 558 G4int nelm = currentMaterial->GetNumberOfElements();
555 559
556 // geometry 560 // geometry
557 G4double sint, cost, phi; 561 G4double sint, cost, phi;
558 G4ThreeVector temp(0.0,0.0,1.0); 562 G4ThreeVector temp(0.0,0.0,1.0);
559 563
560 // current position and direction relative t 564 // current position and direction relative to the end point
561 // because of magnetic field geometry is com 565 // because of magnetic field geometry is computed relatively to the
562 // end point of the step 566 // end point of the step
563 G4ThreeVector dir(0.0,0.0,1.0); 567 G4ThreeVector dir(0.0,0.0,1.0);
564 fDisplacement.set(0.0,0.0,-zPathLength); 568 fDisplacement.set(0.0,0.0,-zPathLength);
565 569
566 G4double mscfac = zPathLength/tPathLength; 570 G4double mscfac = zPathLength/tPathLength;
567 571
568 // start a loop 572 // start a loop
569 G4double x2 = x0; 573 G4double x2 = x0;
570 G4double step, z; 574 G4double step, z;
571 G4bool singleScat; 575 G4bool singleScat;
572 /* 576 /*
573 G4cout << "Start of the loop x1(mm)= " << 577 G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
574 << " 1-cost1= " << 1 - cosThetaMin << " SS 578 << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode
575 << " xtsec= " << xtsec << " Nst= " << 579 << " xtsec= " << xtsec << " Nst= " << nMscSteps << G4endl;
576 */ 580 */
577 do { 581 do {
578 582
579 //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= 583 //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
580 // single scattering case 584 // single scattering case
581 if(singleScatteringMode && x1 > x2) { 585 if(singleScatteringMode && x1 > x2) {
582 fDisplacement += x2*mscfac*dir; 586 fDisplacement += x2*mscfac*dir;
583 break; 587 break;
584 } 588 }
585 589
586 // what is next single of multiple? 590 // what is next single of multiple?
587 if(x1 <= x2) { 591 if(x1 <= x2) {
588 step = x1; 592 step = x1;
589 singleScat = true; 593 singleScat = true;
590 } else { 594 } else {
591 step = x2; 595 step = x2;
592 singleScat = false; 596 singleScat = false;
593 } 597 }
594 598
595 //G4cout << "# step(mm)= "<< step<< " sin 599 //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
596 600
597 // new position 601 // new position
598 fDisplacement += step*mscfac*dir; 602 fDisplacement += step*mscfac*dir;
599 603
600 if(singleScat) { 604 if(singleScat) {
601 605
602 // select element 606 // select element
603 std::size_t i = 0; << 607 G4int i = 0;
604 if(nelm > 1) { 608 if(nelm > 1) {
605 G4double qsec = rndmEngine->flat()*xts << 609 G4double qsec = G4UniformRand()*xtsec;
606 for (; i<nelm; ++i) { if(xsecn[i] >= q << 610 for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
607 } 611 }
608 G4double cosTetM = 612 G4double cosTetM =
609 wokvi->SetupTarget((*theElementVector) << 613 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
610 //G4cout << "!!! " << cosThetaMin << " 614 //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
611 // << prob[i] << G4endl; 615 // << prob[i] << G4endl;
612 temp = wokvi->SampleSingleScattering(cos 616 temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
613 617
614 // direction is changed 618 // direction is changed
615 temp.rotateUz(dir); 619 temp.rotateUz(dir);
616 dir = temp; 620 dir = temp;
617 //G4cout << dir << G4endl; 621 //G4cout << dir << G4endl;
618 622
619 // new proposed step length 623 // new proposed step length
620 x2 -= step; 624 x2 -= step;
621 x1 = -G4Log(rndmEngine->flat())/xtsec; << 625 x1 = -G4Log(G4UniformRand())/xtsec;
622 626
623 // multiple scattering 627 // multiple scattering
624 } else { 628 } else {
625 --nMscSteps; 629 --nMscSteps;
626 x1 -= step; 630 x1 -= step;
627 x2 = x0; 631 x2 = x0;
628 632
629 // sample z in interval 0 - 1 633 // sample z in interval 0 - 1
630 G4bool isFirst = true; 634 G4bool isFirst = true;
631 if(prob2 > 0.0 && rndmEngine->flat() < p << 635 if(prob2 > 0.0 && G4UniformRand() < prob2) { isFirst = false; }
632 do { 636 do {
633 if(isFirst) { z = -G4Log(rndmEngine->f << 637 //z = -z0*G4Log(1.0 - (1.0 - zzz)*G4UniformRand());
634 else { z = G4RandGamma::shoot(r << 638 if(isFirst) { z = -G4Log(G4UniformRand()); }
635 z *= z0; << 639 else { z = G4RandGamma::shoot(2.0,2.0); }
636 // Loop checking, 03-Aug-2015, Vladimi << 640 z *= z0;
637 } while(z > 1.0); 641 } while(z > 1.0);
638 642
639 cost = 1.0 - 2.0*z/*factCM*/; 643 cost = 1.0 - 2.0*z/*factCM*/;
640 if(cost > 1.0) { cost = 1.0; } 644 if(cost > 1.0) { cost = 1.0; }
641 else if(cost < -1.0) { cost =-1.0; } 645 else if(cost < -1.0) { cost =-1.0; }
642 sint = sqrt((1.0 - cost)*(1.0 + cost)); 646 sint = sqrt((1.0 - cost)*(1.0 + cost));
643 phi = twopi*rndmEngine->flat(); << 647 phi = twopi*G4UniformRand();
644 G4double vx1 = sint*std::cos(phi); << 648 G4double vx1 = sint*cos(phi);
645 G4double vy1 = sint*std::sin(phi); << 649 G4double vy1 = sint*sin(phi);
646 650
647 // lateral displacement 651 // lateral displacement
648 if (latDisplasment) { 652 if (latDisplasment) {
649 G4double rms = z0 > 0.0 ? invsqrt12*st << 653 G4double rms = invsqrt12*sqrt(2*z0);
650 G4double r = x0*mscfac; << 654 G4double r = x0*mscfac;
651 G4double dx = r*(0.5*vx1 + rms*G4RandG << 655 G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
652 G4double dy = r*(0.5*vy1 + rms*G4RandG << 656 G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
653 G4double d = r*r - dx*dx - dy*dy; << 657 G4double dz;
654 << 658 G4double d = r*r - dx*dx - dy*dy;
655 // change position << 659
656 if(d >= 0.0) { << 660 // change position
657 temp.set(dx, dy, std::sqrt(d) - r); << 661 if(d >= 0.0) {
658 temp.rotateUz(dir); << 662 dz = sqrt(d) - r;
659 fDisplacement += temp; << 663 temp.set(dx,dy,dz);
660 } << 664 temp.rotateUz(dir);
>> 665 fDisplacement += temp;
>> 666 }
661 } 667 }
662 // change direction 668 // change direction
663 temp.set(vx1,vy1,cost); 669 temp.set(vx1,vy1,cost);
664 temp.rotateUz(dir); 670 temp.rotateUz(dir);
665 dir = temp; 671 dir = temp;
666 } 672 }
667 // Loop checking, 03-Aug-2015, Vladimir Iv <<
668 } while (0 < nMscSteps); 673 } while (0 < nMscSteps);
669 674
670 dir.rotateUz(oldDirection); 675 dir.rotateUz(oldDirection);
671 676
672 //G4cout<<"G4WentzelVIModel sampling is done 677 //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
673 // end of sampling ------------------------- 678 // end of sampling -------------------------------
674 679
675 fParticleChange->ProposeMomentumDirection(di 680 fParticleChange->ProposeMomentumDirection(dir);
676 681
677 // lateral displacement 682 // lateral displacement
678 fDisplacement.rotateUz(oldDirection); 683 fDisplacement.rotateUz(oldDirection);
679 684
680 /* 685 /*
681 G4cout << " r(mm)= " << fDisplacement << 686 G4cout << " r(mm)= " << fDisplacement.mag()
682 << " safety= " << safety << 687 << " safety= " << safety
683 << " trueStep(mm)= " << tPathL << 688 << " trueStep(mm)= " << tPathLength
684 << " geomStep(mm)= " << zPathL << 689 << " geomStep(mm)= " << zPathLength
685 << " x= " << fDisplacement.x() << 690 << " x= " << fDisplacement.x()
686 << " y= " << fDisplacement.y() << 691 << " y= " << fDisplacement.y()
687 << " z= " << fDisplacement.z() << 692 << " z= " << fDisplacement.z()
688 << G4endl; << 693 << G4endl;
689 */ 694 */
690 695
691 //G4cout<< "G4WentzelVIModel::SampleScatteri 696 //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
692 return fDisplacement; 697 return fDisplacement;
693 } 698 }
694 699
695 //....oooOO0OOooo........oooOO0OOooo........oo 700 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
696 701
697 G4double G4WentzelVIModel::ComputeTransportXSe 702 G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume(G4double cosTheta)
698 { 703 {
699 // prepare recomputation of x-sections 704 // prepare recomputation of x-sections
700 const G4ElementVector* theElementVector = cu 705 const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
701 const G4double* theAtomNumDensityVector = 706 const G4double* theAtomNumDensityVector =
702 currentMaterial->GetVecNbOfAtomsPerVolume( 707 currentMaterial->GetVecNbOfAtomsPerVolume();
703 G4int nelm = (G4int)currentMaterial->GetNumb << 708 G4int nelm = currentMaterial->GetNumberOfElements();
704 if(nelm > nelments) { 709 if(nelm > nelments) {
705 nelments = nelm; 710 nelments = nelm;
706 xsecn.resize(nelm); 711 xsecn.resize(nelm);
707 prob.resize(nelm); 712 prob.resize(nelm);
708 } 713 }
709 714
710 // check consistency 715 // check consistency
711 xtsec = 0.0; 716 xtsec = 0.0;
712 if(cosTetMaxNuc >= cosTheta) { return 0.0; } 717 if(cosTetMaxNuc >= cosTheta) { return 0.0; }
713 718
714 G4double cut = (*currentCuts)[currentMateria 719 G4double cut = (*currentCuts)[currentMaterialIndex];
715 if(fixedCut > 0.0) { cut = fixedCut; } 720 if(fixedCut > 0.0) { cut = fixedCut; }
716 721
717 // loop over elements 722 // loop over elements
718 G4double xs = 0.0; 723 G4double xs = 0.0;
719 for (G4int i=0; i<nelm; ++i) { 724 for (G4int i=0; i<nelm; ++i) {
720 G4double costm = 725 G4double costm =
721 wokvi->SetupTarget((*theElementVector)[i << 726 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
722 G4double density = theAtomNumDensityVector 727 G4double density = theAtomNumDensityVector[i];
723 728
724 G4double esec = 0.0; 729 G4double esec = 0.0;
725 if(costm < cosTheta) { 730 if(costm < cosTheta) {
726 731
727 // recompute the transport x-section 732 // recompute the transport x-section
728 if(1.0 > cosTheta) { 733 if(1.0 > cosTheta) {
729 xs += density*wokvi->ComputeTransportC << 734 xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta);
730 } 735 }
731 // recompute the total x-section 736 // recompute the total x-section
732 G4double nucsec = wokvi->ComputeNuclearC 737 G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm);
733 esec = wokvi->ComputeElectronCrossSectio 738 esec = wokvi->ComputeElectronCrossSection(cosTheta, costm);
734 nucsec += esec; 739 nucsec += esec;
735 if(nucsec > 0.0) { esec /= nucsec; } 740 if(nucsec > 0.0) { esec /= nucsec; }
736 xtsec += nucsec*density; 741 xtsec += nucsec*density;
737 } 742 }
738 xsecn[i] = xtsec; 743 xsecn[i] = xtsec;
739 prob[i] = esec; 744 prob[i] = esec;
740 //G4cout << i << " xs= " << xs << " xtsec 745 //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
741 // << " 1-cosTheta= " << 1-cosTheta 746 // << " 1-cosTheta= " << 1-cosTheta
742 // << " 1-cosTetMaxNuc2= " <<1-c << 747 // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
743 } 748 }
744 749
745 //G4cout << "ComputeXS result: xsec(1/mm)= 750 //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
746 // << " txsec(1/mm)= " << xtsec <<G4 << 751 // << " txsec(1/mm)= " << xtsec <<G4endl;
747 return xs; 752 return xs;
748 } 753 }
749 754
750 //....oooOO0OOooo........oooOO0OOooo........oo 755 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
751 756
752 G4double G4WentzelVIModel::ComputeSecondMoment << 757 G4double G4WentzelVIModel:: ComputeSecondMoment(const G4ParticleDefinition* p,
753 G4double kinEnergy) << 758 G4double kinEnergy)
754 { 759 {
755 G4double xs = 0.0; 760 G4double xs = 0.0;
756 761
757 SetupParticle(p); 762 SetupParticle(p);
758 cosTetMaxNuc = wokvi->SetupKinematic(kinEner 763 cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
759 764
760 if(cosTetMaxNuc >= 1.0) { return xs; } 765 if(cosTetMaxNuc >= 1.0) { return xs; }
761 766
762 const G4ElementVector* theElementVector = cu 767 const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
763 const G4double* theAtomNumDensityVector = 768 const G4double* theAtomNumDensityVector =
764 currentMaterial->GetVecNbOfAtomsPerVolume( 769 currentMaterial->GetVecNbOfAtomsPerVolume();
765 std::size_t nelm = currentMaterial->GetNumbe << 770 G4int nelm = currentMaterial->GetNumberOfElements();
766 771
767 G4double cut = (*currentCuts)[currentMateria 772 G4double cut = (*currentCuts)[currentMaterialIndex];
768 if(fixedCut > 0.0) { cut = fixedCut; } 773 if(fixedCut > 0.0) { cut = fixedCut; }
769 774
770 // loop over elements 775 // loop over elements
771 for (std::size_t i=0; i<nelm; ++i) { << 776 for (G4int i=0; i<nelm; ++i) {
772 G4double costm = 777 G4double costm =
773 wokvi->SetupTarget((*theElementVector)[i << 778 wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
774 xs += theAtomNumDensityVector[i] 779 xs += theAtomNumDensityVector[i]
775 *wokvi->ComputeSecondTransportMoment(cos 780 *wokvi->ComputeSecondTransportMoment(costm);
776 } 781 }
777 return xs; 782 return xs;
778 } 783 }
779 784
780 //....oooOO0OOooo........oooOO0OOooo........oo 785 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
781 786
782 void G4WentzelVIModel::SetSingleScatteringFact 787 void G4WentzelVIModel::SetSingleScatteringFactor(G4double val)
783 { 788 {
784 if(val > 0.05) { 789 if(val > 0.05) {
785 ssFactor = val; 790 ssFactor = val;
786 invssFactor = 1.0/(val - 0.05); 791 invssFactor = 1.0/(val - 0.05);
787 } 792 }
788 } 793 }
789 794
790 //....oooOO0OOooo........oooOO0OOooo........oo 795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
791 796