00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104 #include "G4UrbanMscModel93.hh"
00105 #include "G4PhysicalConstants.hh"
00106 #include "G4SystemOfUnits.hh"
00107 #include "Randomize.hh"
00108 #include "G4Electron.hh"
00109 #include "G4LossTableManager.hh"
00110 #include "G4ParticleChangeForMSC.hh"
00111
00112 #include "G4Poisson.hh"
00113 #include "globals.hh"
00114
00115
00116
00117 using namespace std;
00118
00119 G4UrbanMscModel93::G4UrbanMscModel93(const G4String& nam)
00120 : G4VMscModel(nam)
00121 {
00122 masslimite = 0.6*MeV;
00123 lambdalimit = 1.*mm;
00124 fr = 0.02;
00125 taubig = 8.0;
00126 tausmall = 1.e-16;
00127 taulim = 1.e-6;
00128 currentTau = taulim;
00129 tlimitminfix = 1.e-6*mm;
00130 stepmin = tlimitminfix;
00131 smallstep = 1.e10;
00132 currentRange = 0. ;
00133 rangeinit = 0.;
00134 tlimit = 1.e10*mm;
00135 tlimitmin = 10.*tlimitminfix;
00136 tgeom = 1.e50*mm;
00137 geombig = 1.e50*mm;
00138 geommin = 1.e-3*mm;
00139 geomlimit = geombig;
00140 presafety = 0.*mm;
00141
00142 y = 0.;
00143
00144 Zold = 0.;
00145 Zeff = 1.;
00146 Z2 = 1.;
00147 Z23 = 1.;
00148 lnZ = 0.;
00149 coeffth1 = 0.;
00150 coeffth2 = 0.;
00151 coeffc1 = 0.;
00152 coeffc2 = 0.;
00153 scr1ini = fine_structure_const*fine_structure_const*
00154 electron_mass_c2*electron_mass_c2/(0.885*0.885*4.*pi);
00155 scr2ini = 3.76*fine_structure_const*fine_structure_const;
00156 scr1 = 0.;
00157 scr2 = 0.;
00158
00159 theta0max = pi/6.;
00160 rellossmax = 0.50;
00161 third = 1./3.;
00162 particle = 0;
00163 theManager = G4LossTableManager::Instance();
00164 firstStep = true;
00165 inside = false;
00166 insideskin = false;
00167
00168 numlim = 0.01;
00169 xsi = 3.;
00170 ea = exp(-xsi);
00171 eaa = 1.-ea ;
00172
00173 skindepth = skin*stepmin;
00174
00175 mass = proton_mass_c2;
00176 charge = ChargeSquare = 1.0;
00177 currentKinEnergy = currentRadLength = lambda0 = lambdaeff = tPathLength
00178 = zPathLength = par1 = par2 = par3 = 0;
00179
00180 currentMaterialIndex = -1;
00181 fParticleChange = 0;
00182 couple = 0;
00183 SetSampleZ(false);
00184 }
00185
00186
00187
00188 G4UrbanMscModel93::~G4UrbanMscModel93()
00189 {}
00190
00191
00192
00193 void G4UrbanMscModel93::Initialise(const G4ParticleDefinition* p,
00194 const G4DataVector&)
00195 {
00196 skindepth = skin*stepmin;
00197
00198
00199 SetParticle(p);
00200
00201 if(p->GetPDGMass() > MeV) {
00202 G4cout << "### WARNING: G4UrbanMscModel93 model is used for "
00203 << p->GetParticleName() << " !!! " << G4endl;
00204 G4cout << "### This model should be used only for e+-"
00205 << G4endl;
00206 }
00207
00208 fParticleChange = GetParticleChangeForMSC(p);
00209 }
00210
00211
00212
00213 G4double G4UrbanMscModel93::ComputeCrossSectionPerAtom(
00214 const G4ParticleDefinition* part,
00215 G4double KineticEnergy,
00216 G4double AtomicNumber,G4double,
00217 G4double, G4double)
00218 {
00219 static const G4double sigmafactor =
00220 twopi*classic_electr_radius*classic_electr_radius;
00221 static const G4double epsfactor = 2.*electron_mass_c2*electron_mass_c2*
00222 Bohr_radius*Bohr_radius/(hbarc*hbarc);
00223 static const G4double epsmin = 1.e-4 , epsmax = 1.e10;
00224
00225 static const G4double Zdat[15] = { 4., 6., 13., 20., 26., 29., 32., 38., 47.,
00226 50., 56., 64., 74., 79., 82. };
00227
00228 static const G4double Tdat[22] = { 100*eV, 200*eV, 400*eV, 700*eV,
00229 1*keV, 2*keV, 4*keV, 7*keV,
00230 10*keV, 20*keV, 40*keV, 70*keV,
00231 100*keV, 200*keV, 400*keV, 700*keV,
00232 1*MeV, 2*MeV, 4*MeV, 7*MeV,
00233 10*MeV, 20*MeV};
00234
00235
00236 static const G4double celectron[15][22] =
00237 {{1.125,1.072,1.051,1.047,1.047,1.050,1.052,1.054,
00238 1.054,1.057,1.062,1.069,1.075,1.090,1.105,1.111,
00239 1.112,1.108,1.100,1.093,1.089,1.087 },
00240 {1.408,1.246,1.143,1.096,1.077,1.059,1.053,1.051,
00241 1.052,1.053,1.058,1.065,1.072,1.087,1.101,1.108,
00242 1.109,1.105,1.097,1.090,1.086,1.082 },
00243 {2.833,2.268,1.861,1.612,1.486,1.309,1.204,1.156,
00244 1.136,1.114,1.106,1.106,1.109,1.119,1.129,1.132,
00245 1.131,1.124,1.113,1.104,1.099,1.098 },
00246 {3.879,3.016,2.380,2.007,1.818,1.535,1.340,1.236,
00247 1.190,1.133,1.107,1.099,1.098,1.103,1.110,1.113,
00248 1.112,1.105,1.096,1.089,1.085,1.098 },
00249 {6.937,4.330,2.886,2.256,1.987,1.628,1.395,1.265,
00250 1.203,1.122,1.080,1.065,1.061,1.063,1.070,1.073,
00251 1.073,1.070,1.064,1.059,1.056,1.056 },
00252 {9.616,5.708,3.424,2.551,2.204,1.762,1.485,1.330,
00253 1.256,1.155,1.099,1.077,1.070,1.068,1.072,1.074,
00254 1.074,1.070,1.063,1.059,1.056,1.052 },
00255 {11.72,6.364,3.811,2.806,2.401,1.884,1.564,1.386,
00256 1.300,1.180,1.112,1.082,1.073,1.066,1.068,1.069,
00257 1.068,1.064,1.059,1.054,1.051,1.050 },
00258 {18.08,8.601,4.569,3.183,2.662,2.025,1.646,1.439,
00259 1.339,1.195,1.108,1.068,1.053,1.040,1.039,1.039,
00260 1.039,1.037,1.034,1.031,1.030,1.036 },
00261 {18.22,10.48,5.333,3.713,3.115,2.367,1.898,1.631,
00262 1.498,1.301,1.171,1.105,1.077,1.048,1.036,1.033,
00263 1.031,1.028,1.024,1.022,1.021,1.024 },
00264 {14.14,10.65,5.710,3.929,3.266,2.453,1.951,1.669,
00265 1.528,1.319,1.178,1.106,1.075,1.040,1.027,1.022,
00266 1.020,1.017,1.015,1.013,1.013,1.020 },
00267 {14.11,11.73,6.312,4.240,3.478,2.566,2.022,1.720,
00268 1.569,1.342,1.186,1.102,1.065,1.022,1.003,0.997,
00269 0.995,0.993,0.993,0.993,0.993,1.011 },
00270 {22.76,20.01,8.835,5.287,4.144,2.901,2.219,1.855,
00271 1.677,1.410,1.224,1.121,1.073,1.014,0.986,0.976,
00272 0.974,0.972,0.973,0.974,0.975,0.987 },
00273 {50.77,40.85,14.13,7.184,5.284,3.435,2.520,2.059,
00274 1.837,1.512,1.283,1.153,1.091,1.010,0.969,0.954,
00275 0.950,0.947,0.949,0.952,0.954,0.963 },
00276 {65.87,59.06,15.87,7.570,5.567,3.650,2.682,2.182,
00277 1.939,1.579,1.325,1.178,1.108,1.014,0.965,0.947,
00278 0.941,0.938,0.940,0.944,0.946,0.954 },
00279 {55.60,47.34,15.92,7.810,5.755,3.767,2.760,2.239,
00280 1.985,1.609,1.343,1.188,1.113,1.013,0.960,0.939,
00281 0.933,0.930,0.933,0.936,0.939,0.949 }};
00282
00283 static const G4double cpositron[15][22] = {
00284 {2.589,2.044,1.658,1.446,1.347,1.217,1.144,1.110,
00285 1.097,1.083,1.080,1.086,1.092,1.108,1.123,1.131,
00286 1.131,1.126,1.117,1.108,1.103,1.100 },
00287 {3.904,2.794,2.079,1.710,1.543,1.325,1.202,1.145,
00288 1.122,1.096,1.089,1.092,1.098,1.114,1.130,1.137,
00289 1.138,1.132,1.122,1.113,1.108,1.102 },
00290 {7.970,6.080,4.442,3.398,2.872,2.127,1.672,1.451,
00291 1.357,1.246,1.194,1.179,1.178,1.188,1.201,1.205,
00292 1.203,1.190,1.173,1.159,1.151,1.145 },
00293 {9.714,7.607,5.747,4.493,3.815,2.777,2.079,1.715,
00294 1.553,1.353,1.253,1.219,1.211,1.214,1.225,1.228,
00295 1.225,1.210,1.191,1.175,1.166,1.174 },
00296 {17.97,12.95,8.628,6.065,4.849,3.222,2.275,1.820,
00297 1.624,1.382,1.259,1.214,1.202,1.202,1.214,1.219,
00298 1.217,1.203,1.184,1.169,1.160,1.151 },
00299 {24.83,17.06,10.84,7.355,5.767,3.707,2.546,1.996,
00300 1.759,1.465,1.311,1.252,1.234,1.228,1.238,1.241,
00301 1.237,1.222,1.201,1.184,1.174,1.159 },
00302 {23.26,17.15,11.52,8.049,6.375,4.114,2.792,2.155,
00303 1.880,1.535,1.353,1.281,1.258,1.247,1.254,1.256,
00304 1.252,1.234,1.212,1.194,1.183,1.170 },
00305 {22.33,18.01,12.86,9.212,7.336,4.702,3.117,2.348,
00306 2.015,1.602,1.385,1.297,1.268,1.251,1.256,1.258,
00307 1.254,1.237,1.214,1.195,1.185,1.179 },
00308 {33.91,24.13,15.71,10.80,8.507,5.467,3.692,2.808,
00309 2.407,1.873,1.564,1.425,1.374,1.330,1.324,1.320,
00310 1.312,1.288,1.258,1.235,1.221,1.205 },
00311 {32.14,24.11,16.30,11.40,9.015,5.782,3.868,2.917,
00312 2.490,1.925,1.596,1.447,1.391,1.342,1.332,1.327,
00313 1.320,1.294,1.264,1.240,1.226,1.214 },
00314 {29.51,24.07,17.19,12.28,9.766,6.238,4.112,3.066,
00315 2.602,1.995,1.641,1.477,1.414,1.356,1.342,1.336,
00316 1.328,1.302,1.270,1.245,1.231,1.233 },
00317 {38.19,30.85,21.76,15.35,12.07,7.521,4.812,3.498,
00318 2.926,2.188,1.763,1.563,1.484,1.405,1.382,1.371,
00319 1.361,1.330,1.294,1.267,1.251,1.239 },
00320 {49.71,39.80,27.96,19.63,15.36,9.407,5.863,4.155,
00321 3.417,2.478,1.944,1.692,1.589,1.480,1.441,1.423,
00322 1.409,1.372,1.330,1.298,1.280,1.258 },
00323 {59.25,45.08,30.36,20.83,16.15,9.834,6.166,4.407,
00324 3.641,2.648,2.064,1.779,1.661,1.531,1.482,1.459,
00325 1.442,1.400,1.354,1.319,1.299,1.272 },
00326 {56.38,44.29,30.50,21.18,16.51,10.11,6.354,4.542,
00327 3.752,2.724,2.116,1.817,1.692,1.554,1.499,1.474,
00328 1.456,1.412,1.364,1.328,1.307,1.282 }};
00329
00330
00331 static const G4double Tlim = 10.*MeV;
00332 static const G4double beta2lim = Tlim*(Tlim+2.*electron_mass_c2)/
00333 ((Tlim+electron_mass_c2)*(Tlim+electron_mass_c2));
00334 static const G4double bg2lim = Tlim*(Tlim+2.*electron_mass_c2)/
00335 (electron_mass_c2*electron_mass_c2);
00336
00337 static const G4double sig0[15] = {
00338 0.2672*barn, 0.5922*barn, 2.653*barn, 6.235*barn,
00339 11.69*barn , 13.24*barn , 16.12*barn, 23.00*barn ,
00340 35.13*barn , 39.95*barn , 50.85*barn, 67.19*barn ,
00341 91.15*barn , 104.4*barn , 113.1*barn};
00342
00343 static const G4double hecorr[15] = {
00344 120.70, 117.50, 105.00, 92.92, 79.23, 74.510, 68.29,
00345 57.39, 41.97, 36.14, 24.53, 10.21, -7.855, -16.84,
00346 -22.30};
00347
00348 G4double sigma;
00349 SetParticle(part);
00350
00351 Z23 = pow(AtomicNumber,2./3.);
00352
00353
00354
00355
00356
00357 G4double eKineticEnergy = KineticEnergy;
00358
00359 if(mass > electron_mass_c2)
00360 {
00361 G4double TAU = KineticEnergy/mass ;
00362 G4double c = mass*TAU*(TAU+2.)/(electron_mass_c2*(TAU+1.)) ;
00363 G4double w = c-2. ;
00364 G4double tau = 0.5*(w+sqrt(w*w+4.*c)) ;
00365 eKineticEnergy = electron_mass_c2*tau ;
00366 }
00367
00368 G4double eTotalEnergy = eKineticEnergy + electron_mass_c2 ;
00369 G4double beta2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
00370 /(eTotalEnergy*eTotalEnergy);
00371 G4double bg2 = eKineticEnergy*(eTotalEnergy+electron_mass_c2)
00372 /(electron_mass_c2*electron_mass_c2);
00373
00374 G4double eps = epsfactor*bg2/Z23;
00375
00376 if (eps<epsmin) sigma = 2.*eps*eps;
00377 else if(eps<epsmax) sigma = log(1.+2.*eps)-2.*eps/(1.+2.*eps);
00378 else sigma = log(2.*eps)-1.+1./eps;
00379
00380 sigma *= ChargeSquare*AtomicNumber*AtomicNumber/(beta2*bg2);
00381
00382
00383 G4double c1,c2,cc1,cc2,corr;
00384
00385
00386 G4int iZ = 14;
00387 while ((iZ>=0)&&(Zdat[iZ]>=AtomicNumber)) iZ -= 1;
00388 if (iZ==14) iZ = 13;
00389 if (iZ==-1) iZ = 0 ;
00390
00391 G4double ZZ1 = Zdat[iZ];
00392 G4double ZZ2 = Zdat[iZ+1];
00393 G4double ratZ = (AtomicNumber-ZZ1)*(AtomicNumber+ZZ1)/
00394 ((ZZ2-ZZ1)*(ZZ2+ZZ1));
00395
00396 if(eKineticEnergy <= Tlim)
00397 {
00398
00399 G4int iT = 21;
00400 while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy)) iT -= 1;
00401 if(iT==21) iT = 20;
00402 if(iT==-1) iT = 0 ;
00403
00404
00405 G4double T = Tdat[iT], E = T + electron_mass_c2;
00406 G4double b2small = T*(E+electron_mass_c2)/(E*E);
00407
00408 T = Tdat[iT+1]; E = T + electron_mass_c2;
00409 G4double b2big = T*(E+electron_mass_c2)/(E*E);
00410 G4double ratb2 = (beta2-b2small)/(b2big-b2small);
00411
00412 if (charge < 0.)
00413 {
00414 c1 = celectron[iZ][iT];
00415 c2 = celectron[iZ+1][iT];
00416 cc1 = c1+ratZ*(c2-c1);
00417
00418 c1 = celectron[iZ][iT+1];
00419 c2 = celectron[iZ+1][iT+1];
00420 cc2 = c1+ratZ*(c2-c1);
00421
00422 corr = cc1+ratb2*(cc2-cc1);
00423
00424 sigma *= sigmafactor/corr;
00425 }
00426 else
00427 {
00428 c1 = cpositron[iZ][iT];
00429 c2 = cpositron[iZ+1][iT];
00430 cc1 = c1+ratZ*(c2-c1);
00431
00432 c1 = cpositron[iZ][iT+1];
00433 c2 = cpositron[iZ+1][iT+1];
00434 cc2 = c1+ratZ*(c2-c1);
00435
00436 corr = cc1+ratb2*(cc2-cc1);
00437
00438 sigma *= sigmafactor/corr;
00439 }
00440 }
00441 else
00442 {
00443 c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2-beta2lim))/bg2;
00444 c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(beta2-beta2lim))/bg2;
00445 if((AtomicNumber >= ZZ1) && (AtomicNumber <= ZZ2))
00446 sigma = c1+ratZ*(c2-c1) ;
00447 else if(AtomicNumber < ZZ1)
00448 sigma = AtomicNumber*AtomicNumber*c1/(ZZ1*ZZ1);
00449 else if(AtomicNumber > ZZ2)
00450 sigma = AtomicNumber*AtomicNumber*c2/(ZZ2*ZZ2);
00451 }
00452 return sigma;
00453
00454 }
00455
00456
00457
00458 void G4UrbanMscModel93::StartTracking(G4Track* track)
00459 {
00460 SetParticle(track->GetDynamicParticle()->GetDefinition());
00461 firstStep = true;
00462 inside = false;
00463 insideskin = false;
00464 tlimit = geombig;
00465 stepmin = tlimitminfix ;
00466 tlimitmin = 10.*stepmin ;
00467 }
00468
00469
00470
00471 G4double G4UrbanMscModel93::ComputeTruePathLengthLimit(
00472 const G4Track& track,
00473 G4double& currentMinimalStep)
00474 {
00475 tPathLength = currentMinimalStep;
00476 const G4DynamicParticle* dp = track.GetDynamicParticle();
00477 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
00478 G4StepStatus stepStatus = sp->GetStepStatus();
00479 couple = track.GetMaterialCutsCouple();
00480 SetCurrentCouple(couple);
00481 currentMaterialIndex = couple->GetIndex();
00482 currentKinEnergy = dp->GetKineticEnergy();
00483 currentRange = GetRange(particle,currentKinEnergy,couple);
00484 lambda0 = GetTransportMeanFreePath(particle,currentKinEnergy);
00485
00486
00487 if(inside || tPathLength < tlimitminfix) {
00488 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00489 }
00490
00491 if(tPathLength > currentRange) { tPathLength = currentRange; }
00492
00493 presafety = sp->GetSafety();
00494
00495
00496
00497
00498
00499
00500
00501 if(currentRange < presafety)
00502 {
00503 inside = true;
00504 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00505 }
00506
00507
00508
00509 if (steppingAlgorithm == fUseDistanceToBoundary)
00510 {
00511
00512 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
00513
00514
00515 if(currentRange < presafety)
00516 {
00517 inside = true;
00518 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00519 }
00520
00521 smallstep += 1.;
00522 insideskin = false;
00523
00524 if(firstStep || stepStatus == fGeomBoundary)
00525 {
00526 rangeinit = currentRange;
00527 if(firstStep) smallstep = 1.e10;
00528 else smallstep = 1.;
00529
00530
00531
00532 G4double rat = currentKinEnergy/MeV ;
00533 rat = 1.e-3/(rat*(10.+rat)) ;
00534
00535 stepmin = rat*lambda0;
00536 skindepth = skin*stepmin;
00537
00538 tlimitmin = 10.*stepmin;
00539 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00540
00541
00542
00543 if((geomlimit < geombig) && (geomlimit > geommin))
00544 {
00545
00546
00547 if((1.-geomlimit/lambda0) > 0.)
00548 geomlimit = -lambda0*log(1.-geomlimit/lambda0)+tlimitmin ;
00549
00550 if(stepStatus == fGeomBoundary)
00551 tgeom = geomlimit/facgeom;
00552 else
00553 tgeom = 2.*geomlimit/facgeom;
00554 }
00555 else
00556 tgeom = geombig;
00557 }
00558
00559
00560
00561 tlimit = facrange*rangeinit;
00562 if(tlimit < facsafety*presafety)
00563 tlimit = facsafety*presafety;
00564
00565
00566 if(tlimit < tlimitmin) tlimit = tlimitmin;
00567
00568 if(tlimit > tgeom) tlimit = tgeom;
00569
00570
00571
00572
00573
00574 if((tPathLength < tlimit) && (tPathLength < presafety) &&
00575 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
00576 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00577
00578
00579 if(smallstep < skin)
00580 {
00581 tlimit = stepmin;
00582 insideskin = true;
00583 }
00584 else if(geomlimit < geombig)
00585 {
00586 if(geomlimit > skindepth)
00587 {
00588 if(tlimit > geomlimit-0.999*skindepth)
00589 tlimit = geomlimit-0.999*skindepth;
00590 }
00591 else
00592 {
00593 insideskin = true;
00594 if(tlimit > stepmin) tlimit = stepmin;
00595 }
00596 }
00597
00598 if(tlimit < stepmin) tlimit = stepmin;
00599
00600
00601 if(firstStep || ((smallstep == skin) && !insideskin))
00602 {
00603 G4double temptlimit = tlimit;
00604 if(temptlimit > tlimitmin)
00605 {
00606 do {
00607 temptlimit = G4RandGauss::shoot(tlimit,0.3*tlimit);
00608 } while ((temptlimit < tlimitmin) ||
00609 (temptlimit > 2.*tlimit-tlimitmin));
00610 }
00611 else
00612 temptlimit = tlimitmin;
00613 if(tPathLength > temptlimit) tPathLength = temptlimit;
00614 }
00615 else
00616 {
00617 if(tPathLength > tlimit) tPathLength = tlimit ;
00618 }
00619
00620 }
00621
00622
00623 else if(steppingAlgorithm == fUseSafety)
00624 {
00625
00626
00627 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
00628 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
00629
00630
00631 if(currentRange < presafety)
00632 {
00633 inside = true;
00634 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00635 }
00636
00637 if(firstStep || stepStatus == fGeomBoundary)
00638 {
00639 rangeinit = currentRange;
00640 fr = facrange;
00641
00642 if(mass < masslimite)
00643 {
00644 if(lambda0 > currentRange)
00645 rangeinit = lambda0;
00646 if(lambda0 > lambdalimit)
00647 fr *= 0.75+0.25*lambda0/lambdalimit;
00648 }
00649
00650
00651 G4double rat = currentKinEnergy/MeV ;
00652 rat = 1.e-3/(rat*(10.+rat)) ;
00653 tlimitmin = 10.*lambda0*rat;
00654 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
00655 }
00656
00657 tlimit = fr*rangeinit;
00658
00659 if(tlimit < facsafety*presafety)
00660 tlimit = facsafety*presafety;
00661
00662
00663 if(tlimit < tlimitmin) tlimit = tlimitmin;
00664
00665 if(tPathLength > tlimit) tPathLength = tlimit;
00666
00667 }
00668
00669
00670 else
00671 {
00672 if (stepStatus == fGeomBoundary)
00673 {
00674 if (currentRange > lambda0) tlimit = facrange*currentRange;
00675 else tlimit = facrange*lambda0;
00676
00677 if(tlimit < tlimitmin) tlimit = tlimitmin;
00678 if(tPathLength > tlimit) tPathLength = tlimit;
00679 }
00680 }
00681
00682
00683 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
00684 }
00685
00686
00687
00688 G4double G4UrbanMscModel93::ComputeGeomPathLength(G4double)
00689 {
00690 firstStep = false;
00691 lambdaeff = lambda0;
00692 par1 = -1. ;
00693 par2 = par3 = 0. ;
00694
00695
00696 zPathLength = tPathLength;
00697
00698
00699 if(tPathLength < tlimitminfix) return zPathLength;
00700
00701
00702
00703 if(tPathLength > currentRange)
00704 tPathLength = currentRange ;
00705
00706 G4double tau = tPathLength/lambda0 ;
00707
00708 if ((tau <= tausmall) || insideskin) {
00709 zPathLength = tPathLength;
00710 if(zPathLength > lambda0) zPathLength = lambda0;
00711 return zPathLength;
00712 }
00713
00714 G4double zmean = tPathLength;
00715 if (tPathLength < currentRange*dtrl) {
00716 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
00717 else zmean = lambda0*(1.-exp(-tau));
00718 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
00719 par1 = 1./currentRange ;
00720 par2 = 1./(par1*lambda0) ;
00721 par3 = 1.+par2 ;
00722 if(tPathLength < currentRange)
00723 zmean = (1.-exp(par3*log(1.-tPathLength/currentRange)))/(par1*par3) ;
00724 else
00725 zmean = 1./(par1*par3) ;
00726 } else {
00727 G4double T1 = GetEnergy(particle,currentRange-tPathLength,couple);
00728 G4double lambda1 = GetTransportMeanFreePath(particle,T1);
00729
00730 par1 = (lambda0-lambda1)/(lambda0*tPathLength) ;
00731 par2 = 1./(par1*lambda0) ;
00732 par3 = 1.+par2 ;
00733 zmean = (1.-exp(par3*log(lambda1/lambda0)))/(par1*par3) ;
00734 }
00735
00736 zPathLength = zmean ;
00737
00738
00739 if(samplez)
00740 {
00741 const G4double ztmax = 0.99 ;
00742 G4double zt = zmean/tPathLength ;
00743
00744 if (tPathLength > stepmin && zt < ztmax)
00745 {
00746 G4double u,cz1;
00747 if(zt >= third)
00748 {
00749 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
00750 cz1 = 1.+cz ;
00751 G4double u0 = cz/cz1 ;
00752 G4double grej ;
00753 do {
00754 u = exp(log(G4UniformRand())/cz1) ;
00755 grej = exp(cz*log(u/u0))*(1.-u)/(1.-u0) ;
00756 } while (grej < G4UniformRand()) ;
00757 }
00758 else
00759 {
00760 cz1 = 1./zt-1.;
00761 u = 1.-exp(log(G4UniformRand())/cz1) ;
00762 }
00763 zPathLength = tPathLength*u ;
00764 }
00765 }
00766
00767 if(zPathLength > lambda0) zPathLength = lambda0;
00768
00769 return zPathLength;
00770 }
00771
00772
00773
00774 G4double G4UrbanMscModel93::ComputeTrueStepLength(G4double geomStepLength)
00775 {
00776
00777 if(geomStepLength == zPathLength && tPathLength <= currentRange)
00778 return tPathLength;
00779
00780
00781 zPathLength = geomStepLength;
00782 tPathLength = geomStepLength;
00783 if(geomStepLength < tlimitminfix) return tPathLength;
00784
00785
00786 if((geomStepLength > lambda0*tausmall) && !insideskin)
00787 {
00788 if(par1 < 0.)
00789 tPathLength = -lambda0*log(1.-geomStepLength/lambda0) ;
00790 else
00791 {
00792 if(par1*par3*geomStepLength < 1.)
00793 tPathLength = (1.-exp(log(1.-par1*par3*geomStepLength)/par3))/par1 ;
00794 else
00795 tPathLength = currentRange;
00796 }
00797 }
00798 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
00799
00800
00801
00802 return tPathLength;
00803 }
00804
00805
00806
00807 G4double G4UrbanMscModel93::ComputeTheta0(G4double trueStepLength,
00808 G4double KineticEnergy)
00809 {
00810
00811
00812
00813 static const G4double c_highland = 13.6*MeV ;
00814 G4double betacp = sqrt(currentKinEnergy*(currentKinEnergy+2.*mass)*
00815 KineticEnergy*(KineticEnergy+2.*mass)/
00816 ((currentKinEnergy+mass)*(KineticEnergy+mass)));
00817 y = trueStepLength/currentRadLength;
00818 G4double theta0 = c_highland*std::abs(charge)*sqrt(y)/betacp;
00819 y = log(y);
00820
00821 G4double corr = coeffth1+coeffth2*y;
00822
00823 theta0 *= corr ;
00824
00825 return theta0;
00826 }
00827
00828
00829
00830 G4ThreeVector&
00831 G4UrbanMscModel93::SampleScattering(const G4ThreeVector& oldDirection,
00832 G4double safety)
00833 {
00834 fDisplacement.set(0.0,0.0,0.0);
00835 G4double kineticEnergy = currentKinEnergy;
00836 if (tPathLength > currentRange*dtrl) {
00837 kineticEnergy = GetEnergy(particle,currentRange-tPathLength,couple);
00838 } else {
00839 kineticEnergy -= tPathLength*GetDEDX(particle,currentKinEnergy,couple);
00840 }
00841 if((kineticEnergy <= eV) || (tPathLength <= tlimitminfix) ||
00842 (tPathLength/tausmall < lambda0)) { return fDisplacement; }
00843
00844 G4double cth = SampleCosineTheta(tPathLength,kineticEnergy);
00845
00846
00847 if(std::fabs(cth) > 1.) { return fDisplacement; }
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867
00868
00869
00870
00871
00872 G4double sth = sqrt((1.0 - cth)*(1.0 + cth));
00873 G4double phi = twopi*G4UniformRand();
00874 G4double dirx = sth*cos(phi);
00875 G4double diry = sth*sin(phi);
00876
00877 G4ThreeVector newDirection(dirx,diry,cth);
00878 newDirection.rotateUz(oldDirection);
00879 fParticleChange->ProposeMomentumDirection(newDirection);
00880
00881 if (latDisplasment && safety > tlimitminfix) {
00882
00883 G4double r = SampleDisplacement();
00884
00885
00886
00887
00888
00889
00890
00891 if(r > 0.)
00892 {
00893 G4double latcorr = LatCorrelation();
00894 if(latcorr > r) latcorr = r;
00895
00896
00897
00898 G4double Phi = 0.;
00899 if(std::abs(r*sth) < latcorr)
00900 Phi = twopi*G4UniformRand();
00901 else
00902 {
00903 G4double psi = std::acos(latcorr/(r*sth));
00904 if(G4UniformRand() < 0.5)
00905 Phi = phi+psi;
00906 else
00907 Phi = phi-psi;
00908 }
00909
00910 dirx = r*std::cos(Phi);
00911 diry = r*std::sin(Phi);
00912
00913 fDisplacement.set(dirx,diry,0.0);
00914 fDisplacement.rotateUz(oldDirection);
00915 }
00916 }
00917 return fDisplacement;
00918 }
00919
00920
00921
00922 G4double G4UrbanMscModel93::SampleCosineTheta(G4double trueStepLength,
00923 G4double KineticEnergy)
00924 {
00925 G4double cth = 1. ;
00926 G4double tau = trueStepLength/lambda0 ;
00927
00928 Zeff = couple->GetMaterial()->GetTotNbOfElectPerVolume()/
00929 couple->GetMaterial()->GetTotNbOfAtomsPerVolume() ;
00930
00931 if(Zold != Zeff)
00932 UpdateCache();
00933
00934 if(insideskin)
00935 {
00936
00937 G4double mean = trueStepLength/stepmin ;
00938
00939 G4int n = G4Poisson(mean);
00940 if(n > 0)
00941 {
00942
00943 G4double mom2 = KineticEnergy*(2.*mass+KineticEnergy);
00944 G4double beta2 = mom2/((KineticEnergy+mass)*(KineticEnergy+mass));
00945 G4double ascr = scr1/mom2;
00946 ascr *= 1.13+scr2/beta2;
00947 G4double ascr1 = 1.+2.*ascr;
00948 G4double bp1=ascr1+1.;
00949 G4double bm1=ascr1-1.;
00950
00951
00952 G4double ct,st,phi;
00953 G4double sx=0.,sy=0.,sz=0.;
00954 for(G4int i=1; i<=n; i++)
00955 {
00956 ct = ascr1-bp1*bm1/(2.*G4UniformRand()+bm1);
00957 if(ct < -1.) ct = -1.;
00958 if(ct > 1.) ct = 1.;
00959 st = sqrt(1.-ct*ct);
00960 phi = twopi*G4UniformRand();
00961 sx += st*cos(phi);
00962 sy += st*sin(phi);
00963 sz += ct;
00964 }
00965 cth = sz/sqrt(sx*sx+sy*sy+sz*sz);
00966 }
00967 }
00968 else
00969 {
00970 if(trueStepLength >= currentRange*dtrl)
00971 {
00972 if(par1*trueStepLength < 1.)
00973 tau = -par2*log(1.-par1*trueStepLength) ;
00974
00975
00976 else if(1.-KineticEnergy/currentKinEnergy > taulim)
00977 tau = taubig ;
00978 }
00979 currentTau = tau ;
00980 lambdaeff = trueStepLength/currentTau;
00981 currentRadLength = couple->GetMaterial()->GetRadlen();
00982
00983 if (tau >= taubig) cth = -1.+2.*G4UniformRand();
00984 else if (tau >= tausmall)
00985 {
00986 G4double xmeanth, x2meanth;
00987 if(tau < numlim) {
00988 xmeanth = 1.0 - tau*(1.0 - 0.5*tau);
00989 x2meanth= 1.0 - tau*(5.0 - 6.25*tau)/3.;
00990 } else {
00991 xmeanth = exp(-tau);
00992 x2meanth = (1.+2.*exp(-2.5*tau))/3.;
00993 }
00994 G4double relloss = 1.-KineticEnergy/currentKinEnergy;
00995
00996 if(relloss > rellossmax)
00997 return SimpleScattering(xmeanth,x2meanth);
00998
00999 G4double theta0 = ComputeTheta0(trueStepLength,KineticEnergy);
01000
01001
01002
01003
01004
01005 G4double theta2 = theta0*theta0;
01006
01007 if(theta2 < tausmall) { return cth; }
01008
01009 if(theta0 > theta0max) {
01010 return SimpleScattering(xmeanth,x2meanth);
01011 }
01012
01013 G4double x = theta2*(1.0 - theta2/12.);
01014 if(theta2 > numlim) {
01015 G4double sth = 2.*sin(0.5*theta0);
01016 x = sth*sth;
01017 }
01018
01019 G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/eaa;
01020 G4double x0 = 1. - xsi*x;
01021
01022
01023
01024 if(xmean1 <= 0.999*xmeanth) {
01025 return SimpleScattering(xmeanth,x2meanth);
01026 }
01027
01028 G4double c = coeffc1+coeffc2*y;
01029
01030
01031 if(c < 1.9) {
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 c = 1.9;
01043 }
01044
01045 if(fabs(c-3.) < 0.001) { c = 3.001; }
01046 else if(fabs(c-2.) < 0.001) { c = 2.001; }
01047
01048 G4double c1 = c-1.;
01049
01050
01051 G4double b = 1.+(c-xsi)*x;
01052
01053 G4double b1 = b+1.;
01054 G4double bx = c*x;
01055
01056 G4double eb1 = pow(b1,c1);
01057 G4double ebx = pow(bx,c1);
01058 G4double d = ebx/eb1;
01059
01060
01061 G4double xmean2 = (x0 + d - (bx - b1*d)/(c-2.))/(1. - d);
01062
01063 G4double f1x0 = ea/eaa;
01064 G4double f2x0 = c1/(c*(1. - d));
01065 G4double prob = f2x0/(f1x0+f2x0);
01066
01067 G4double qprob = xmeanth/(prob*xmean1+(1.-prob)*xmean2);
01068
01069
01070
01071
01072
01073 if(G4UniformRand() < qprob)
01074 {
01075 G4double var = 0;
01076 if(G4UniformRand() < prob) {
01077 cth = 1.+log(ea+G4UniformRand()*eaa)*x;
01078 } else {
01079 var = (1.0 - d)*G4UniformRand();
01080 if(var < numlim*d) {
01081 var /= (d*c1);
01082 cth = -1.0 + var*(1.0 - 0.5*var*c)*(2. + (c - xsi)*x);
01083 } else {
01084 cth = 1. + x*(c - xsi - c*pow(var + d, -1.0/c1));
01085
01086 }
01087 }
01088 if(KineticEnergy > 5*GeV && cth < 0.9) {
01089 G4cout << "G4UrbanMscModel93::SampleCosineTheta: E(GeV)= "
01090 << KineticEnergy/GeV
01091 << " 1-cosT= " << 1 - cth
01092 << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
01093 << " tau= " << tau
01094 << " prob= " << prob << " var= " << var << G4endl;
01095 G4cout << " c= " << c << " qprob= " << qprob << " eb1= " << eb1
01096 << " ebx= " << ebx
01097 << " c1= " << c1 << " b= " << b << " b1= " << b1
01098 << " bx= " << bx << " d= " << d
01099 << " ea= " << ea << " eaa= " << eaa << G4endl;
01100 }
01101 }
01102 else {
01103 cth = -1.+2.*G4UniformRand();
01104 if(KineticEnergy > 5*GeV) {
01105 G4cout << "G4UrbanMscModel93::SampleCosineTheta: E(GeV)= "
01106 << KineticEnergy/GeV
01107 << " length(mm)= " << trueStepLength << " Zeff= " << Zeff
01108 << " qprob= " << qprob << G4endl;
01109 }
01110 }
01111 }
01112 }
01113 return cth ;
01114 }
01115
01116
01117
01118 G4double G4UrbanMscModel93::SimpleScattering(G4double xmeanth,G4double x2meanth)
01119 {
01120
01121
01122 G4double a = (2.*xmeanth+9.*x2meanth-3.)/(2.*xmeanth-3.*x2meanth+1.);
01123 G4double prob = (a+2.)*xmeanth/a;
01124
01125
01126 G4double cth = 1.;
01127 if(G4UniformRand() < prob)
01128 cth = -1.+2.*exp(log(G4UniformRand())/(a+1.));
01129 else
01130 cth = -1.+2.*G4UniformRand();
01131 return cth;
01132 }
01133
01134
01135
01136 G4double G4UrbanMscModel93::SampleDisplacement()
01137 {
01138
01139 const G4double kappa = 2.5;
01140 const G4double kappapl1 = kappa+1.;
01141 const G4double kappami1 = kappa-1.;
01142
01143 G4double rmean = 0.0;
01144 if ((currentTau >= tausmall) && !insideskin) {
01145 if (currentTau < taulim) {
01146 rmean = kappa*currentTau*currentTau*currentTau*
01147 (1.-kappapl1*currentTau*0.25)/6. ;
01148
01149 } else {
01150 G4double etau = 0.0;
01151 if (currentTau<taubig) etau = exp(-currentTau);
01152 rmean = -kappa*currentTau;
01153 rmean = -exp(rmean)/(kappa*kappami1);
01154 rmean += currentTau-kappapl1/kappa+kappa*etau/kappami1;
01155 }
01156 if (rmean>0.) rmean = 2.*lambdaeff*sqrt(rmean/3.0);
01157 else rmean = 0.;
01158 }
01159
01160 if(rmean == 0.) return rmean;
01161
01162
01163 G4double rmax = (tPathLength-zPathLength)*(tPathLength+zPathLength);
01164 if(rmax <= 0.)
01165 rmax = 0.;
01166 else
01167 rmax = sqrt(rmax);
01168
01169 if(rmean >= rmax) return rmax;
01170
01171 return rmean;
01172
01173
01174
01175
01176
01177
01178
01179
01180
01181
01182
01183
01184
01185
01186 }
01187
01188
01189
01190 G4double G4UrbanMscModel93::LatCorrelation()
01191 {
01192 const G4double kappa = 2.5;
01193 const G4double kappami1 = kappa-1.;
01194
01195 G4double latcorr = 0.;
01196 if((currentTau >= tausmall) && !insideskin)
01197 {
01198 if(currentTau < taulim)
01199 latcorr = lambdaeff*kappa*currentTau*currentTau*
01200 (1.-(kappa+1.)*currentTau/3.)/3.;
01201 else
01202 {
01203 G4double etau = 0.;
01204 if(currentTau < taubig) etau = exp(-currentTau);
01205 latcorr = -kappa*currentTau;
01206 latcorr = exp(latcorr)/kappami1;
01207 latcorr += 1.-kappa*etau/kappami1 ;
01208 latcorr *= 2.*lambdaeff/3. ;
01209 }
01210 }
01211
01212 return latcorr;
01213 }
01214
01215