Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
G4UrbanMscModel96.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class file
30//
31//
32// File name: G4UrbanMscModel96
33//
34// Author: Laszlo Urban
35//
36// Creation date: 26.09.2012
37//
38// Created from G4UrbanMscModel95
39//
40// New parametrization for theta0
41// Correction for very small step length
42//
43// Class Description:
44//
45// Implementation of the model of multiple scattering based on
46// H.W.Lewis Phys Rev 78 (1950) 526 and others
47
48// -------------------------------------------------------------------
49// In its present form the model can be used for simulation
50// of the e-/e+ multiple scattering
51//
52
53
54//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57#include "G4UrbanMscModel96.hh"
59#include "G4SystemOfUnits.hh"
60#include "Randomize.hh"
61#include "G4Electron.hh"
62#include "G4LossTableManager.hh"
64
65#include "G4Poisson.hh"
66#include "globals.hh"
67
68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69
70using namespace std;
71
73 : G4VMscModel(nam)
74{
75 masslimite = 0.6*MeV;
76 lambdalimit = 1.*mm;
77 fr = 0.02;
78 taubig = 8.0;
79 tausmall = 1.e-16;
80 taulim = 1.e-6;
81 currentTau = taulim;
82 tlimitminfix = 1.e-6*mm;
83 stepmin = tlimitminfix;
84 smallstep = 1.e10;
85 currentRange = 0. ;
86 rangeinit = 0.;
87 tlimit = 1.e10*mm;
88 tlimitmin = 10.*tlimitminfix;
89 tgeom = 1.e50*mm;
90 geombig = 1.e50*mm;
91 geommin = 1.e-3*mm;
92 geomlimit = geombig;
93 presafety = 0.*mm;
94 //facsafety = 0.50 ;
95
96 y = 0.;
97 z = 0.;
98
99 Zold = 0.;
100 Zeff = 1.;
101 Z2 = 1.;
102 Z23 = 1.;
103 lnZ = 0.;
104
105 coeffc1 = 0.;
106 coeffc2 = 0.;
107 coeffc3 = 0.;
108 coeffc4 = 0.;
109 scr1ini = fine_structure_const*fine_structure_const*
110 electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
111 scr2ini = 3.76*fine_structure_const*fine_structure_const;
112 scr1 = 0.;
113 scr2 = 0.;
114
115 theta0max = pi/6.;
116 rellossmax = 0.50;
117 third = 1./3.;
118 particle = 0;
119 theManager = G4LossTableManager::Instance();
120 firstStep = true;
121 inside = false;
122 insideskin = false;
123
124 skindepth = skin*stepmin;
125
126 mass = proton_mass_c2;
127 charge = ChargeSquare = 1.0;
128 currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
129 = zPathLength = par1 = par2 = par3 = 0;
130
131 currentMaterialIndex = -1;
132 fParticleChange = 0;
133 couple = 0;
134 SetSampleZ(true);
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
138
140{}
141
142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
143
145 const G4DataVector&)
146{
147 skindepth = skin*stepmin;
148 // trackID = -1;
149
150 // set values of some data members
151 SetParticle(p);
152
153 if(p->GetPDGMass() > MeV) {
154 G4cout << "### WARNING: G4UrbanMscModel96 model is used for "
155 << p->GetParticleName() << " !!! " << G4endl;
156 G4cout << "### This model should be used only for e+-"
157 << G4endl;
158 }
159
160 fParticleChange = GetParticleChangeForMSC(p);
161
162 //samplez = true;
163 //G4cout << "### G4UrbanMscModel96::Initialise done!" << G4endl;
164}
165
166//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
167
169 const G4ParticleDefinition* part,
170 G4double KineticEnergy,
171 G4double AtomicNumber,G4double,
173{
174 const G4double sigmafactor = twopi*classic_electr_radius*classic_electr_radius;
175 const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
176 Bohr_radius*Bohr_radius/(hbarc*hbarc);
177 const G4double epsmin = 1.e-4 , epsmax = 1.e10;
178
179 const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38., 47.,
180 50., 56., 64., 74., 79., 82. };
181
182 const G4double Tdat[22] = { 100*eV, 200*eV, 400*eV, 700*eV,
183 1*keV, 2*keV, 4*keV, 7*keV,
184 10*keV, 20*keV, 40*keV, 70*keV,
185 100*keV, 200*keV, 400*keV, 700*keV,
186 1*MeV, 2*MeV, 4*MeV, 7*MeV,
187 10*MeV, 20*MeV};
188
189 // corr. factors for e-/e+ lambda for T <= Tlim
190 G4double celectron[15][22] =
191 {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
192 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
193 1.112,1.108,1.100,1.093,1.089,1.087 },
194 {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
195 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
196 1.109,1.105,1.097,1.090,1.086,1.082 },
197 {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
198 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
199 1.131,1.124,1.113,1.104,1.099,1.098 },
200 {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
201 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
202 1.112,1.105,1.096,1.089,1.085,1.098 },
203 {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
204 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
205 1.073,1.070,1.064,1.059,1.056,1.056 },
206 {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
207 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
208 1.074,1.070,1.063,1.059,1.056,1.052 },
209 {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
210 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
211 1.068,1.064,1.059,1.054,1.051,1.050 },
212 {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
213 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
214 1.039,1.037,1.034,1.031,1.030,1.036 },
215 {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
216 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
217 1.031,1.028,1.024,1.022,1.021,1.024 },
218 {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
219 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
220 1.020,1.017,1.015,1.013,1.013,1.020 },
221 {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
222 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
223 0.995,0.993,0.993,0.993,0.993,1.011 },
224 {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
225 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
226 0.974,0.972,0.973,0.974,0.975,0.987 },
227 {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
228 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
229 0.950,0.947,0.949,0.952,0.954,0.963 },
230 {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
231 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
232 0.941,0.938,0.940,0.944,0.946,0.954 },
233 {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
234 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
235 0.933,0.930,0.933,0.936,0.939,0.949 }};
236
237 G4double cpositron[15][22] = {
238 {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
239 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
240 1.131,1.126,1.117,1.108,1.103,1.100 },
241 {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
242 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
243 1.138,1.132,1.122,1.113,1.108,1.102 },
244 {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
245 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
246 1.203,1.190,1.173,1.159,1.151,1.145 },
247 {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
248 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
249 1.225,1.210,1.191,1.175,1.166,1.174 },
250 {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
251 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
252 1.217,1.203,1.184,1.169,1.160,1.151 },
253 {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
254 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
255 1.237,1.222,1.201,1.184,1.174,1.159 },
256 {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
257 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
258 1.252,1.234,1.212,1.194,1.183,1.170 },
259 {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
260 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
261 1.254,1.237,1.214,1.195,1.185,1.179 },
262 {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
263 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
264 1.312,1.288,1.258,1.235,1.221,1.205 },
265 {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
266 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
267 1.320,1.294,1.264,1.240,1.226,1.214 },
268 {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
269 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
270 1.328,1.302,1.270,1.245,1.231,1.233 },
271 {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
272 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
273 1.361,1.330,1.294,1.267,1.251,1.239 },
274 {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
275 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
276 1.409,1.372,1.330,1.298,1.280,1.258 },
277 {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
278 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
279 1.442,1.400,1.354,1.319,1.299,1.272 },
280 {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
281 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
282 1.456,1.412,1.364,1.328,1.307,1.282 }};
283
284 //data/corrections for T > Tlim
285 G4double Tlim = 10.*MeV;
286 G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
287 ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
288 G4double bg2lim = Tlim*(Tlim+2.*electron_mass_c2)/
289 (electron_mass_c2*electron_mass_c2);
290
291 G4double sig0[15] = {0.2672*barn, 0.5922*barn, 2.653*barn, 6.235*barn,
292 11.69*barn , 13.24*barn , 16.12*barn, 23.00*barn ,
293 35.13*barn , 39.95*barn , 50.85*barn, 67.19*barn ,
294 91.15*barn , 104.4*barn , 113.1*barn};
295
296 G4double hecorr[15] = {120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
297 57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
298 -22.30};
299
300 G4double sigma;
301 SetParticle(part);
302
303 Z23 = pow(AtomicNumber,2./3.);
304
305 // correction if particle .ne. e-/e+
306 // compute equivalent kinetic energy
307 // lambda depends on p*beta ....
308
309 G4double eKineticEnergy = KineticEnergy;
310
311 if(mass > electron_mass_c2)
312 {
313 G4double TAU = KineticEnergy/mass ;
314 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
315 G4double w = c-2. ;
316 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
317 eKineticEnergy = electron_mass_c2*tau ;
318 }
319
320 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
321 G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
322 /(eTotalEnergy*eTotalEnergy);
323 G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
324 /(electron_mass_c2*electron_mass_c2);
325
326 G4double eps = epsfactor*bg2/Z23;
327
328 if (eps<epsmin) sigma = 2.*eps*eps;
329 else if(eps<epsmax) sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
330 else sigma = log(2.*eps)-1.+1./eps;
331
332 sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
333
334 // interpolate in AtomicNumber and beta2
335 G4double c1,c2,cc1,cc2,corr;
336
337 // get bin number in Z
338 G4int iZ = 14;
339 while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
340 if (iZ==14) iZ = 13;
341 if (iZ==-1) iZ = 0 ;
342
343 G4double ZZ1 = Zdat[iZ];
344 G4double ZZ2 = Zdat[iZ+1];
345 G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
346 ((ZZ2-ZZ1)*(ZZ2+ZZ1));
347
348 if(eKineticEnergy <= Tlim)
349 {
350 // get bin number in T (beta2)
351 G4int iT = 21;
352 while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
353 if(iT==21) iT = 20;
354 if(iT==-1) iT = 0 ;
355
356 // calculate betasquare values
357 G4double T = Tdat[iT], E = T + electron_mass_c2;
358 G4double b2small = T*(E+electron_mass_c2)/(E*E);
359
360 T = Tdat[iT+1]; E = T + electron_mass_c2;
361 G4double b2big = T*(E+electron_mass_c2)/(E*E);
362 G4double ratb2 = (beta2-b2small)/(b2big-b2small);
363
364 if (charge < 0.)
365 {
366 c1 = celectron[iZ][iT];
367 c2 = celectron[iZ+1][iT];
368 cc1 = c1+ratZ*(c2-c1);
369
370 c1 = celectron[iZ][iT+1];
371 c2 = celectron[iZ+1][iT+1];
372 cc2 = c1+ratZ*(c2-c1);
373
374 corr = cc1+ratb2*(cc2-cc1);
375
376 sigma *= sigmafactor/corr;
377 }
378 else
379 {
380 c1 = cpositron[iZ][iT];
381 c2 = cpositron[iZ+1][iT];
382 cc1 = c1+ratZ*(c2-c1);
383
384 c1 = cpositron[iZ][iT+1];
385 c2 = cpositron[iZ+1][iT+1];
386 cc2 = c1+ratZ*(c2-c1);
387
388 corr = cc1+ratb2*(cc2-cc1);
389
390 sigma *= sigmafactor/corr;
391 }
392 }
393 else
394 {
395 c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
396 c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
397 if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
398 sigma = c1+ratZ*(c2-c1) ;
399 else if(AtomicNumber < ZZ1)
400 sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
401 else if(AtomicNumber > ZZ2)
402 sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
403 }
404 return sigma;
405
406}
407
408//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
409
411{
412 SetParticle(track->GetDynamicParticle()->GetDefinition());
413 firstStep = true;
414 inside = false;
415 insideskin = false;
416 tlimit = geombig;
417}
418
419//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
420
422 const G4Track& track,
423 G4double& currentMinimalStep)
424{
425 tPathLength = currentMinimalStep;
426 const G4DynamicParticle* dp = track.GetDynamicParticle();
427
428 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
429 G4StepStatus stepStatus = sp->GetStepStatus();
430 couple = track.GetMaterialCutsCouple();
431 SetCurrentCouple(couple);
432 currentMaterialIndex = couple->GetIndex();
433 currentKinEnergy = dp->GetKineticEnergy();
434 currentRange = GetRange(particle,currentKinEnergy,couple);
435 lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
436
437 if(tPathLength > currentRange) { tPathLength = currentRange; }
438
439 // stop here if small range particle
440 if(inside || tPathLength < tlimitminfix) {
441 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
442 }
443
444 presafety = sp->GetSafety();
445 /*
446 G4cout << "G4Urban96::StepLimit tPathLength= "
447 <<tPathLength<<" safety= " << presafety
448 << " range= " <<currentRange<< " lambda= "<<lambda0
449 << " Alg: " << steppingAlgorithm <<G4endl;
450 */
451 // far from geometry boundary
452 if(currentRange < presafety)
453 {
454 inside = true;
455 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
456 }
457
458 // standard version
459 //
461 {
462 //compute geomlimit and presafety
463 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
464
465 // is it far from boundary ?
466 if(currentRange < presafety)
467 {
468 inside = true;
469 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
470 }
471
472 smallstep += 1.;
473 insideskin = false;
474
475 if(firstStep || (stepStatus == fGeomBoundary))
476 {
477 rangeinit = currentRange;
478 if(firstStep) smallstep = 1.e10;
479 else smallstep = 1.;
480
481 //define stepmin here (it depends on lambda!)
482 //rough estimation of lambda_elastic/lambda_transport
483 G4double rat = currentKinEnergy/MeV ;
484 rat = 1.e-3/(rat*(10.+rat)) ;
485 //stepmin ~ lambda_elastic
486 stepmin = rat*lambda0;
487 skindepth = skin*stepmin;
488 //define tlimitmin
489 tlimitmin = 10.*stepmin;
490 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
491 //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
492 // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
493 // constraint from the geometry
494 if((geomlimit < geombig) && (geomlimit > geommin))
495 {
496 // geomlimit is a geometrical step length
497 // transform it to true path length (estimation)
498 if((1.-geomlimit/lambda0) > 0.)
499 geomlimit = -lambda0*log(1.-geomlimit/lambda0)+tlimitmin ;
500
501 if(stepStatus == fGeomBoundary)
502 tgeom = geomlimit/facgeom;
503 else
504 tgeom = 2.*geomlimit/facgeom;
505 }
506 else
507 tgeom = geombig;
508 }
509
510
511 //step limit
512 tlimit = facrange*rangeinit;
513
514 //lower limit for tlimit
515 if(tlimit < tlimitmin) tlimit = tlimitmin;
516
517 if(tlimit > tgeom) tlimit = tgeom;
518
519 //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
520 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
521
522 // shortcut
523 if((tPathLength < tlimit) && (tPathLength < presafety) &&
524 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
525 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
526
527 // step reduction near to boundary
528 if(smallstep < skin)
529 {
530 tlimit = stepmin;
531 insideskin = true;
532 }
533 else if(geomlimit < geombig)
534 {
535 if(geomlimit > skindepth)
536 {
537 if(tlimit > geomlimit-0.999*skindepth)
538 tlimit = geomlimit-0.999*skindepth;
539 }
540 else
541 {
542 insideskin = true;
543 if(tlimit > stepmin) tlimit = stepmin;
544 }
545 }
546
547 if(tlimit < stepmin) tlimit = stepmin;
548
549 // randomize 1st step or 1st 'normal' step in volume
550 if(firstStep || ((smallstep == skin) && !insideskin))
551 {
552 G4double temptlimit = tlimit;
553 if(temptlimit > tlimitmin)
554 {
555 do {
556 temptlimit = G4RandGauss::shoot(tlimit,0.3*tlimit);
557 } while ((temptlimit < tlimitmin) ||
558 (temptlimit > 2.*tlimit-tlimitmin));
559 }
560 else
561 temptlimit = tlimitmin;
562 if(tPathLength > temptlimit) tPathLength = temptlimit;
563 }
564 else
565 {
566 if(tPathLength > tlimit) tPathLength = tlimit ;
567 }
568
569 }
570 // for 'normal' simulation with or without magnetic field
571 // there no small step/single scattering at boundaries
572 else if(steppingAlgorithm == fUseSafety)
573 {
574 // compute presafety again if presafety <= 0 and no boundary
575 // i.e. when it is needed for optimization purposes
576 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
577 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
578 /*
579 G4cout << "presafety= " << presafety
580 << " firstStep= " << firstStep
581 << " stepStatus= " << stepStatus
582 << G4endl;
583 */
584 // is far from boundary
585 if(currentRange < presafety)
586 {
587 inside = true;
588 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
589 }
590
591 if(firstStep || stepStatus == fGeomBoundary)
592 {
593 rangeinit = currentRange;
594 fr = facrange;
595 // 9.1 like stepping for e+/e- only (not for muons,hadrons)
596 if(mass < masslimite)
597 {
598 if(lambda0 > currentRange)
599 rangeinit = lambda0;
600 if(lambda0 > lambdalimit)
601 fr *= 0.75+0.25*lambda0/lambdalimit;
602 }
603
604 //lower limit for tlimit
605 G4double rat = currentKinEnergy/MeV ;
606 rat = 1.e-3/(rat*(10.+rat)) ;
607 tlimitmin = 10.*lambda0*rat;
608 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
609 }
610 //step limit
611 tlimit = fr*rangeinit;
612
613 if(tlimit < facsafety*presafety)
614 tlimit = facsafety*presafety;
615
616 //lower limit for tlimit
617 if(tlimit < tlimitmin) tlimit = tlimitmin;
618
619 if(tPathLength > tlimit) tPathLength = tlimit;
620
621 }
622
623 // version similar to 7.1 (needed for some experiments)
624 else
625 {
626 if (stepStatus == fGeomBoundary)
627 {
628 if (currentRange > lambda0) tlimit = facrange*currentRange;
629 else tlimit = facrange*lambda0;
630
631 if(tlimit < tlimitmin) tlimit = tlimitmin;
632 if(tPathLength > tlimit) tPathLength = tlimit;
633 }
634 }
635 //G4cout << "tPathLength= " << tPathLength
636 // << " currentMinimalStep= " << currentMinimalStep << G4endl;
637 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
638}
639
640//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
641
643{
644 firstStep = false;
645 lambdaeff = lambda0;
646 par1 = -1. ;
647 par2 = par3 = 0. ;
648
649 // do the true -> geom transformation
650 zPathLength = tPathLength;
651
652 // z = t for very small tPathLength
653 if(tPathLength < tlimitminfix) return zPathLength;
654
655 // this correction needed to run MSC with eIoni and eBrem inactivated
656 // and makes no harm for a normal run
657 if(tPathLength > currentRange)
658 tPathLength = currentRange ;
659
660 G4double tau = tPathLength/lambda0 ;
661
662 if ((tau <= tausmall) || insideskin) {
663 zPathLength = tPathLength;
664 if(zPathLength > lambda0) zPathLength = lambda0;
665 return zPathLength;
666 }
667
668 G4double zmean = tPathLength;
669 if (tPathLength < currentRange*dtrl) {
670 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
671 else zmean = lambda0*(1.-exp(-tau));
672 zPathLength = zmean ;
673 return zPathLength;
674
675 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
676 par1 = 1./currentRange ;
677 par2 = 1./(par1*lambda0) ;
678 par3 = 1.+par2 ;
679 if(tPathLength < currentRange)
680 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
681 else {
682 zmean = 1./(par1*par3) ;
683 }
684 zPathLength = zmean ;
685 return zPathLength;
686
687 } else {
688 G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
689 G4double lambda1 = GetTransportMeanFreePath(particle,T1);
690
691 par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
692 par2 = 1./(par1*lambda0) ;
693 par3 = 1.+par2 ;
694 zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
695 }
696
697 zPathLength = zmean ;
698
699 // sample z
700 if(samplez)
701 {
702 const G4double ztmax = 0.999 ;
703 G4double zt = zmean/tPathLength ;
704
705 if (tPathLength > stepmin && zt < ztmax)
706 {
707 G4double u,cz1;
708 if(zt >= third)
709 {
710 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
711 cz1 = 1.+cz ;
712 G4double u0 = cz/cz1 ;
713 G4double grej ;
714 do {
715 u = exp(log(G4UniformRand())/cz1) ;
716 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
717 } while (grej < G4UniformRand()) ;
718 }
719 else
720 {
721 u = 2.*zt*G4UniformRand();
722 }
723 if(u > 1.0) { u = 1.0; }
724 zPathLength = tPathLength*u ;
725 }
726 }
727
728 if(zPathLength > lambda0) { zPathLength = lambda0; }
729 //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda0 << G4endl;
730 return zPathLength;
731}
732
733//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
734
736{
737 // step defined other than transportation
738 if(geomStepLength == zPathLength)
739 { return tPathLength; }
740
741 zPathLength = geomStepLength;
742
743 // t = z for very small step
744 if(geomStepLength < tlimitminfix) {
745 tPathLength = geomStepLength;
746
747 // recalculation
748 } else {
749
750 G4double tlength = geomStepLength;
751 if((geomStepLength > lambda0*tausmall) && !insideskin) {
752
753 if(par1 < 0.) {
754 tlength = -lambda0*log(1.-geomStepLength/lambda0) ;
755 } else {
756 if(par1*par3*geomStepLength < 1.) {
757 tlength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
758 } else {
759 tlength = currentRange;
760 }
761 }
762 if(tlength < geomStepLength) { tlength = geomStepLength; }
763 else if(tlength > tPathLength) { tlength = tPathLength; }
764 }
765 tPathLength = tlength;
766 }
767 //G4cout << "Urban96::ComputeTrueLength: tPathLength= " << tPathLength
768 // << " step= " << geomStepLength << G4endl;
769 return tPathLength;
770}
771
772//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
773
775 G4double KineticEnergy)
776{
777 // for all particles take the width of the central part
778 // from a parametrization similar to the Highland formula
779 // ( Highland formula: Particle Physics Booklet, July 2002, eq. 26.10)
780 const G4double c_highland = 13.6*MeV ;
781 G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
782 KineticEnergy*(KineticEnergy+2.*mass)/
783 ((currentKinEnergy+mass)*(KineticEnergy+mass)));
784 y = trueStepLength/currentRadLength;
785 G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
786 y = log(y);
787 z = log(currentTau)-y;
788
789 // correction factor from e- scattering data
790 G4double corr = 1.0925*exp(8.05922e-2*y)*(1.+1.06221e-2*z);
791 if(y > -5.5) corr += 3.47658e-2*exp(3.*log(y+5.5))*(1.+1.33868*z);
792
793
794 theta0 *= corr ;
795
796 return theta0;
797}
798
799//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
800
803 G4double safety)
804{
805 fDisplacement.set(0.0,0.0,0.0);
806 G4double kineticEnergy = currentKinEnergy;
807 if (tPathLength > currentRange*dtrl) {
808 kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
809 } else {
810 kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
811 }
812
813 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix) ||
814 (tPathLength/tausmall < lambda0)) { return fDisplacement; }
815
816 G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
817
818 // protection against 'bad' cth values
819 if(std::fabs(cth) > 1.) { return fDisplacement; }
820
821 // extra protection agaist high energy particles backscattered
822 if(cth < 1.0 - 1000*tPathLength/lambda0 && kineticEnergy > 20*MeV) {
823 //G4cout << "Warning: large scattering E(MeV)= " << kineticEnergy
824 // << " s(mm)= " << tPathLength/mm
825 // << " 1-cosTheta= " << 1.0 - cth << G4endl;
826 // do Gaussian central scattering
827 if(kineticEnergy > GeV && cth < 0.0) {
829 ed << dynParticle->GetDefinition()->GetParticleName()
830 << " E(MeV)= " << kineticEnergy/MeV
831 << " Step(mm)= " << tPathLength/mm
832 << " in " << CurrentCouple()->GetMaterial()->GetName()
833 << " CosTheta= " << cth
834 << " is too big - the angle is resampled" << G4endl;
835 G4Exception("G4UrbanMscModel96::SampleScattering","em0004",
836 JustWarning, ed,"");
837 }
838 do {
839 cth = 1.0 + 2*log(G4UniformRand())*tPathLength/lambda0;
840 } while(cth < -1.0);
841 }
842
843 G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
844 G4double phi = twopi*G4UniformRand();
845 G4double dirx = sth*cos(phi);
846 G4double diry = sth*sin(phi);
847
848 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
849 G4ThreeVector newDirection(dirx,diry,cth);
850 newDirection.rotateUz(oldDirection);
851 fParticleChange->ProposeMomentumDirection(newDirection);
852 /*
853 G4cout << "G4UrbanMscModel96::SampleSecondaries: e(MeV)= " << kineticEnergy
854 << " sinTheta= " << sth << " safety(mm)= " << safety
855 << " trueStep(mm)= " << tPathLength
856 << " geomStep(mm)= " << zPathLength
857 << G4endl;
858 */
859 if (latDisplasment && safety > tlimitminfix) {
860
861 G4double r = SampleDisplacement();
862 /*
863 G4cout << "G4UrbanMscModel96::SampleSecondaries: e(MeV)= " << kineticEnergy
864 << " sinTheta= " << sth << " r(mm)= " << r
865 << " trueStep(mm)= " << tPathLength
866 << " geomStep(mm)= " << zPathLength
867 << G4endl;
868 */
869 if(r > 0.)
870 {
871 G4double latcorr = LatCorrelation();
872 if(latcorr > r) latcorr = r;
873
874 // sample direction of lateral displacement
875 // compute it from the lateral correlation
876 G4double Phi = 0.;
877 if(std::abs(r*sth) < latcorr)
878 Phi = twopi*G4UniformRand();
879 else
880 {
881 G4double psi = std::acos(latcorr/(r*sth));
882 if(G4UniformRand() < 0.5)
883 Phi = phi+psi;
884 else
885 Phi = phi-psi;
886 }
887
888 dirx = std::cos(Phi);
889 diry = std::sin(Phi);
890
891 fDisplacement.set(r*dirx,r*diry,0.0);
892 fDisplacement.rotateUz(oldDirection);
893 }
894 }
895 return fDisplacement;
896}
897
898//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
899
900G4double G4UrbanMscModel96::SampleCosineTheta(G4double trueStepLength,
901 G4double KineticEnergy)
902{
903 G4double cth = 1. ;
904 G4double tau = trueStepLength/lambda0;
905 currentTau = tau;
906 lambdaeff = lambda0;
907
908 Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
910
911 if(Zold != Zeff)
912 UpdateCache();
913
914 if(insideskin)
915 {
916 //no scattering, single or plural scattering
917 G4double mean = trueStepLength/stepmin ;
918
919 G4int n = G4Poisson(mean);
920 if(n > 0)
921 {
922 //screening (Moliere-Bethe)
923 G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
924 G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
925 G4double ascr = scr1/mom2;
926 ascr *= 1.13+scr2/beta2;
927 G4double ascr1 = 1.+2.*ascr;
928 G4double bp1=ascr1+1.;
929 G4double bm1=ascr1-1.;
930
931 // single scattering from screened Rutherford x-section
932 G4double ct,st,phi;
933 G4double sx=0.,sy=0.,sz=0.;
934 for(G4int i=1; i<=n; i++)
935 {
936 ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
937 if(ct < -1.) ct = -1.;
938 if(ct > 1.) ct = 1.;
939 st = sqrt(1.-ct*ct);
940 phi = twopi*G4UniformRand();
941 sx += st*cos(phi);
942 sy += st*sin(phi);
943 sz += ct;
944 }
945 cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
946 }
947 }
948 else
949 {
950 G4double lambda1 = GetTransportMeanFreePath(particle,KineticEnergy);
951 if(std::fabs(lambda1/lambda0 - 1) > 0.01 && lambda1 > 0.)
952 {
953 // mean tau value
954 tau = trueStepLength*log(lambda0/lambda1)/(lambda0-lambda1);
955 }
956
957 currentTau = tau ;
958 lambdaeff = trueStepLength/currentTau;
959 currentRadLength = couple->GetMaterial()->GetRadlen();
960
961 if (tau >= taubig) cth = -1.+2.*G4UniformRand();
962 else if (tau >= tausmall)
963 {
964 // correction for very small step length
965 G4double rat = currentKinEnergy/MeV ;
966 rat = 1.e-3/(rat*(10.+rat)) ;
967 G4double lelastic = rat*lambda0;
968 if(trueStepLength < 5.*lelastic)
969 if(G4UniformRand() < exp(-trueStepLength/lelastic))
970 return cth;
971
972 G4double x0 = 1.;
973 G4double a = 1., ea = 0., eaa = 1.;
974 G4double b=1.,b1=2.,bx=1.,eb1=3.,ebx=1.;
975 G4double prob = 1. ;
976 G4double xmean1 = 1., xmean2 = 0.;
977 G4double xmeanth = exp(-tau);
978 G4double x2meanth = (1.+2.*exp(-2.5*tau))/3.;
979
980 G4double relloss = 1.-KineticEnergy/currentKinEnergy;
981 if(relloss > rellossmax)
982 return SimpleScattering(xmeanth,x2meanth);
983
984 G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
985
986 // protection for very small angles
987 if(theta0*theta0 < tausmall) return cth;
988
989 if(theta0 > theta0max)
990 return SimpleScattering(xmeanth,x2meanth);
991 G4double sth = sin(0.5*theta0);
992 a = 0.25/(sth*sth);
993
994 // parameter for tail
995 G4double u = exp(log(tau)/6.);
996 G4double x = log(lambdaeff/currentRadLength);
997 G4double xsi = coeffc1+u*(coeffc2+coeffc3*u)+coeffc4*x;
998 G4double c = xsi;
999
1000 if(abs(c-3.) < 0.001) c = 3.001;
1001 if(abs(c-2.) < 0.001) c = 2.001;
1002 if(abs(c-1.) < 0.001) c = 1.001;
1003
1004 ea = exp(-xsi);
1005 eaa = 1.-ea ;
1006 xmean1 = 1.-(1.-(1.+xsi)*ea)/(a*eaa);
1007 x0 = 1.-xsi/a;
1008
1009 if(xmean1 <= 0.999*xmeanth)
1010 return SimpleScattering(xmeanth,x2meanth);
1011
1012 G4double c1 = c-1.;
1013
1014 //from continuity of derivatives
1015 b = 1.+(c-xsi)/a;
1016
1017 b1 = b+1.;
1018 bx = c/a;
1019 eb1 = exp(c1*log(b1));
1020 ebx = exp(c1*log(bx));
1021
1022 xmean2 = (x0*eb1+ebx-(eb1*bx-b1*ebx)/(c-2.))/(eb1-ebx);
1023
1024 G4double f1x0 = a*ea/eaa;
1025 G4double f2x0 = c1*eb1/(bx*(eb1-ebx));
1026 prob = f2x0/(f1x0+f2x0);
1027
1028 // sampling of costheta
1029 G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
1030 if(G4UniformRand() < qprob)
1031 {
1032 if(G4UniformRand() < prob)
1033 cth = 1.+log(ea+G4UniformRand()*eaa)/a ;
1034 else
1035 cth = b-b1*bx/exp(log(ebx+(eb1-ebx)*G4UniformRand())/c1) ;
1036 }
1037 else
1038 cth = 2.*G4UniformRand()-1.;
1039 }
1040 }
1041 return cth ;
1042}
1043
1044//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1045
1046G4double G4UrbanMscModel96::SimpleScattering(G4double xmeanth,G4double x2meanth)
1047{
1048 // 'large angle scattering'
1049 // 2 model functions with correct xmean and x2mean
1050 G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
1051 G4double prob = (a+2.)*xmeanth/a;
1052
1053 // sampling
1054 G4double cth = 1.;
1055 if(G4UniformRand() < prob)
1056 cth = -1.+2.*exp(log(G4UniformRand())/(a+1.));
1057 else
1058 cth = -1.+2.*G4UniformRand();
1059 return cth;
1060}
1061
1062//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1063
1064G4double G4UrbanMscModel96::SampleDisplacement()
1065{
1066 G4double r = 0.0;
1067 if ((currentTau >= tausmall) && !insideskin) {
1068 G4double rmax = sqrt((tPathLength-zPathLength)*(tPathLength+zPathLength));
1069 r = rmax*exp(log(G4UniformRand())/3.);
1070 }
1071 return r;
1072}
1073
1074//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1075
1076G4double G4UrbanMscModel96::LatCorrelation()
1077{
1078 const G4double kappa = 2.5;
1079 const G4double kappami1 = kappa-1.;
1080
1081 G4double latcorr = 0.;
1082 if((currentTau >= tausmall) && !insideskin)
1083 {
1084 if(currentTau < taulim)
1085 latcorr = lambdaeff*kappa*currentTau*currentTau*
1086 (1.-(kappa+1.)*currentTau/3.)/3.;
1087 else
1088 {
1089 G4double etau = 0.;
1090 if(currentTau < taubig) etau = exp(-currentTau);
1091 latcorr = -kappa*currentTau;
1092 latcorr = exp(latcorr)/kappami1;
1093 latcorr += 1.-kappa*etau/kappami1 ;
1094 latcorr *= 2.*lambdaeff/3. ;
1095 }
1096 }
1097
1098 return latcorr;
1099}
1100
1101//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@ JustWarning
@ fUseSafety
@ fUseDistanceToBoundary
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:50
G4StepStatus
Definition: G4StepStatus.hh:51
@ fGeomBoundary
Definition: G4StepStatus.hh:54
double G4double
Definition: G4Types.hh:64
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:53
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
const G4ThreeVector & GetMomentumDirection() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
static G4LossTableManager * Instance()
const G4Material * GetMaterial() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:208
G4double GetTotNbOfElectPerVolume() const
Definition: G4Material.hh:211
G4double GetRadlen() const
Definition: G4Material.hh:219
const G4String & GetName() const
Definition: G4Material.hh:177
void ProposeMomentumDirection(const G4ThreeVector &Pfinal)
const G4String & GetParticleName() const
G4StepPoint * GetPreStepPoint() const
const G4DynamicParticle * GetDynamicParticle() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
const G4Step * GetStep() const
G4double ComputeGeomPathLength(G4double truePathLength)
G4ThreeVector & SampleScattering(const G4DynamicParticle *, G4double safety)
G4double ComputeTrueStepLength(G4double geomStepLength)
G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)
void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4double ComputeTheta0(G4double truePathLength, G4double KineticEnergy)
G4UrbanMscModel96(const G4String &nam="UrbanMsc96")
G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *particle, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=0., G4double emax=DBL_MAX)
void StartTracking(G4Track *)
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:370
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:377
G4double dtrl
Definition: G4VMscModel.hh:180
G4double GetDEDX(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:273
G4double facrange
Definition: G4VMscModel.hh:176
G4double ComputeGeomLimit(const G4Track &, G4double &presafety, G4double limit)
Definition: G4VMscModel.hh:256
G4bool samplez
Definition: G4VMscModel.hh:188
G4double skin
Definition: G4VMscModel.hh:179
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:332
G4double GetEnergy(const G4ParticleDefinition *part, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:304
G4double GetRange(const G4ParticleDefinition *part, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4VMscModel.hh:288
G4double ComputeSafety(const G4ThreeVector &position, G4double limit)
Definition: G4VMscModel.hh:238
G4MscStepLimitType steppingAlgorithm
Definition: G4VMscModel.hh:186
G4ParticleChangeForMSC * GetParticleChangeForMSC(const G4ParticleDefinition *p=0)
Definition: G4VMscModel.cc:89
G4double ConvertTrueToGeom(G4double &tLength, G4double &gLength)
Definition: G4VMscModel.hh:246
G4bool latDisplasment
Definition: G4VMscModel.hh:189
G4double facsafety
Definition: G4VMscModel.hh:178
G4ThreeVector fDisplacement
Definition: G4VMscModel.hh:185
G4double facgeom
Definition: G4VMscModel.hh:177
void SetSampleZ(G4bool)
Definition: G4VMscModel.hh:231
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76