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