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