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