Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // ------------------------------------------- 27 // 28 // File name: G4UrbanAdjointMscModel 29 // 30 // Author: Laszlo Urban 31 // 32 // Creation date: 19.02.2013 33 // 34 // Created from G4UrbanAdjointMscModel96 35 // 36 // Implementation of the model of multiple sca 37 // H.W.Lewis Phys Rev 78 (1950) 526 and others 38 // In its present form the model can be used 39 // of the e-/e+ multiple scattering 40 // ------------------------------------------- 41 42 #include "G4UrbanAdjointMscModel.hh" 43 44 #include "globals.hh" 45 #include "G4Electron.hh" 46 #include "G4Exp.hh" 47 #include "G4Log.hh" 48 #include "G4LossTableManager.hh" 49 #include "G4ParticleChangeForMSC.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4Poisson.hh" 52 #include "G4Positron.hh" 53 #include "G4Pow.hh" 54 #include "G4SystemOfUnits.hh" 55 #include "Randomize.hh" 56 57 //....oooOO0OOooo........oooOO0OOooo........oo 58 59 using namespace std; 60 61 //....oooOO0OOooo........oooOO0OOooo........oo 62 63 G4UrbanAdjointMscModel::G4UrbanAdjointMscModel 64 : G4VMscModel(nam) 65 { 66 masslimite = 0.6 * MeV; 67 lambdalimit = 1. * mm; 68 fr = 0.02; 69 taubig = 8.0; 70 tausmall = 1.e-16; 71 taulim = 1.e-6; 72 currentTau = taulim; 73 tlimitminfix = 0.01 * nm; 74 tlimitminfix2 = 1. * nm; 75 stepmin = tlimitminfix; 76 smallstep = 1.e10; 77 currentRange = 0.; 78 rangeinit = 0.; 79 tlimit = 1.e10 * mm; 80 tlimitmin = 10. * tlimitminfix; 81 tgeom = 1.e50 * mm; 82 geombig = 1.e50 * mm; 83 geommin = 1.e-3 * mm; 84 geomlimit = geombig; 85 presafety = 0. * mm; 86 87 facsafety = 0.6; 88 89 Zold = 0.; 90 Zeff = 1.; 91 Z2 = 1.; 92 Z23 = 1.; 93 lnZ = 0.; 94 coeffth1 = 0.; 95 coeffth2 = 0.; 96 coeffc1 = 0.; 97 coeffc2 = 0.; 98 coeffc3 = 0.; 99 coeffc4 = 0.; 100 particle = nullptr; 101 102 positron = G4Positron::Positron(); 103 theManager = G4LossTableManager::Instance 104 rndmEngineMod = G4Random::getTheEngine(); 105 106 firstStep = true; 107 insideskin = false; 108 latDisplasmentbackup = false; 109 displacementFlag = true; 110 111 rangecut = geombig; 112 drr = 0.35; 113 finalr = 10. * um; 114 115 skindepth = skin * stepmin; 116 117 mass = proton_mass_c2; 118 charge = ChargeSquare = 1.0; 119 currentKinEnergy = currentRadLength = lambda 120 zPathLength = par1 = par2 = par3 = 0.; 121 122 currentMaterialIndex = -1; 123 fParticleChange = nullptr; 124 couple = nullptr; 125 } 126 127 //....oooOO0OOooo........oooOO0OOooo........oo 128 G4UrbanAdjointMscModel::~G4UrbanAdjointMscMode 129 130 //....oooOO0OOooo........oooOO0OOooo........oo 131 void G4UrbanAdjointMscModel::Initialise(const 132 const 133 { 134 const G4ParticleDefinition* p1 = p; 135 136 if(p->GetParticleName() == "adj_e-") 137 p1 = G4Electron::Electron(); 138 // set values of some data members 139 SetParticle(p1); 140 141 fParticleChange = GetParticleChangeForMSC(p1 142 143 latDisplasmentbackup = latDisplasment; 144 } 145 146 //....oooOO0OOooo........oooOO0OOooo........oo 147 148 G4double G4UrbanAdjointMscModel::ComputeCrossS 149 const G4ParticleDefinition* part, G4double K 150 G4double AtomicNumber, G4double, G4double, G 151 { 152 static constexpr G4double epsmin = 1.e-4; 153 static constexpr G4double epsmax = 1.e10; 154 155 static constexpr G4double Zdat[15] = { 4., 156 47., 157 158 // corr. factors for e-/e+ lambda for T <= T 159 static constexpr G4double celectron[15][22] 160 { 1.125, 1.072, 1.051, 1.047, 1.047, 1.050 161 1.054, 1.057, 1.062, 1.069, 1.075, 1.090 162 1.112, 1.108, 1.100, 1.093, 1.089, 1.087 163 { 1.408, 1.246, 1.143, 1.096, 1.077, 1.059 164 1.052, 1.053, 1.058, 1.065, 1.072, 1.087 165 1.109, 1.105, 1.097, 1.090, 1.086, 1.082 166 { 2.833, 2.268, 1.861, 1.612, 1.486, 1.309 167 1.136, 1.114, 1.106, 1.106, 1.109, 1.119 168 1.131, 1.124, 1.113, 1.104, 1.099, 1.098 169 { 3.879, 3.016, 2.380, 2.007, 1.818, 1.535 170 1.190, 1.133, 1.107, 1.099, 1.098, 1.103 171 1.112, 1.105, 1.096, 1.089, 1.085, 1.098 172 { 6.937, 4.330, 2.886, 2.256, 1.987, 1.628 173 1.203, 1.122, 1.080, 1.065, 1.061, 1.063 174 1.073, 1.070, 1.064, 1.059, 1.056, 1.056 175 { 9.616, 5.708, 3.424, 2.551, 2.204, 1.762 176 1.256, 1.155, 1.099, 1.077, 1.070, 1.068 177 1.074, 1.070, 1.063, 1.059, 1.056, 1.052 178 { 11.72, 6.364, 3.811, 2.806, 2.401, 1.884 179 1.300, 1.180, 1.112, 1.082, 1.073, 1.066 180 1.068, 1.064, 1.059, 1.054, 1.051, 1.050 181 { 18.08, 8.601, 4.569, 3.183, 2.662, 2.025 182 1.339, 1.195, 1.108, 1.068, 1.053, 1.040 183 1.039, 1.037, 1.034, 1.031, 1.030, 1.036 184 { 18.22, 10.48, 5.333, 3.713, 3.115, 2.367 185 1.498, 1.301, 1.171, 1.105, 1.077, 1.048 186 1.031, 1.028, 1.024, 1.022, 1.021, 1.024 187 { 14.14, 10.65, 5.710, 3.929, 3.266, 2.453 188 1.528, 1.319, 1.178, 1.106, 1.075, 1.040 189 1.020, 1.017, 1.015, 1.013, 1.013, 1.020 190 { 14.11, 11.73, 6.312, 4.240, 3.478, 2.566 191 1.569, 1.342, 1.186, 1.102, 1.065, 1.022 192 0.995, 0.993, 0.993, 0.993, 0.993, 1.011 193 { 22.76, 20.01, 8.835, 5.287, 4.144, 2.901 194 1.677, 1.410, 1.224, 1.121, 1.073, 1.014 195 0.974, 0.972, 0.973, 0.974, 0.975, 0.987 196 { 50.77, 40.85, 14.13, 7.184, 5.284, 3.435 197 1.837, 1.512, 1.283, 1.153, 1.091, 1.010 198 0.950, 0.947, 0.949, 0.952, 0.954, 0.963 199 { 65.87, 59.06, 15.87, 7.570, 5.567, 3.650 200 1.939, 1.579, 1.325, 1.178, 1.108, 1.014 201 0.941, 0.938, 0.940, 0.944, 0.946, 0.954 202 { 55.60, 47.34, 15.92, 7.810, 5.755, 3.767 203 1.985, 1.609, 1.343, 1.188, 1.113, 1.013 204 0.933, 0.930, 0.933, 0.936, 0.939, 0.949 205 }; 206 207 static constexpr G4double cpositron[15][22] 208 { 2.589, 2.044, 1.658, 1.446, 1.347, 1.217 209 1.097, 1.083, 1.080, 1.086, 1.092, 1.108 210 1.131, 1.126, 1.117, 1.108, 1.103, 1.100 211 { 3.904, 2.794, 2.079, 1.710, 1.543, 1.325 212 1.122, 1.096, 1.089, 1.092, 1.098, 1.114 213 1.138, 1.132, 1.122, 1.113, 1.108, 1.102 214 { 7.970, 6.080, 4.442, 3.398, 2.872, 2.127 215 1.357, 1.246, 1.194, 1.179, 1.178, 1.188 216 1.203, 1.190, 1.173, 1.159, 1.151, 1.145 217 { 9.714, 7.607, 5.747, 4.493, 3.815, 2.777 218 1.553, 1.353, 1.253, 1.219, 1.211, 1.214 219 1.225, 1.210, 1.191, 1.175, 1.166, 1.174 220 { 17.97, 12.95, 8.628, 6.065, 4.849, 3.222 221 1.624, 1.382, 1.259, 1.214, 1.202, 1.202 222 1.217, 1.203, 1.184, 1.169, 1.160, 1.151 223 { 24.83, 17.06, 10.84, 7.355, 5.767, 3.707 224 1.759, 1.465, 1.311, 1.252, 1.234, 1.228 225 1.237, 1.222, 1.201, 1.184, 1.174, 1.159 226 { 23.26, 17.15, 11.52, 8.049, 6.375, 4.114 227 1.880, 1.535, 1.353, 1.281, 1.258, 1.247 228 1.252, 1.234, 1.212, 1.194, 1.183, 1.170 229 { 22.33, 18.01, 12.86, 9.212, 7.336, 4.702 230 2.015, 1.602, 1.385, 1.297, 1.268, 1.251 231 1.254, 1.237, 1.214, 1.195, 1.185, 1.179 232 { 33.91, 24.13, 15.71, 10.80, 8.507, 5.467 233 2.407, 1.873, 1.564, 1.425, 1.374, 1.330 234 1.312, 1.288, 1.258, 1.235, 1.221, 1.205 235 { 32.14, 24.11, 16.30, 11.40, 9.015, 5.782 236 2.490, 1.925, 1.596, 1.447, 1.391, 1.342 237 1.320, 1.294, 1.264, 1.240, 1.226, 1.214 238 { 29.51, 24.07, 17.19, 12.28, 9.766, 6.238 239 2.602, 1.995, 1.641, 1.477, 1.414, 1.356 240 1.328, 1.302, 1.270, 1.245, 1.231, 1.233 241 { 38.19, 30.85, 21.76, 15.35, 12.07, 7.521 242 2.926, 2.188, 1.763, 1.563, 1.484, 1.405 243 1.361, 1.330, 1.294, 1.267, 1.251, 1.239 244 { 49.71, 39.80, 27.96, 19.63, 15.36, 9.407 245 3.417, 2.478, 1.944, 1.692, 1.589, 1.480 246 1.409, 1.372, 1.330, 1.298, 1.280, 1.258 247 { 59.25, 45.08, 30.36, 20.83, 16.15, 9.834 248 3.641, 2.648, 2.064, 1.779, 1.661, 1.531 249 1.442, 1.400, 1.354, 1.319, 1.299, 1.272 250 { 56.38, 44.29, 30.50, 21.18, 16.51, 10.11 251 3.752, 2.724, 2.116, 1.817, 1.692, 1.554 252 1.456, 1.412, 1.364, 1.328, 1.307, 1.282 253 }; 254 255 // data/corrections for T > Tlim 256 static constexpr G4double hecorr[15] = { 120 257 79. 258 41. 259 -7. 260 261 G4double sigma; 262 SetParticle(part); 263 264 Z23 = G4Pow::GetInstance()->Z23(G4lrint(Atom 265 266 // correction if particle .ne. e-/e+ 267 // compute equivalent kinetic energy 268 // lambda depends on p*beta .... 269 270 G4double eKineticEnergy = KineticEnergy; 271 272 if(mass > electron_mass_c2) 273 { 274 G4double tau1 = KineticEnergy / mass; 275 G4double c = mass * tau1 * (tau1 + 2.) / 276 G4double w = c - 2.; 277 G4double tau = 0.5 * (w + sqrt(w * w + 4. 278 eKineticEnergy = electron_mass_c2 * tau; 279 } 280 281 G4double eTotalEnergy = eKineticEnergy + ele 282 G4double beta2 = eKineticEnergy * (eT 283 (eTotalEnergy * eTotalEnerg 284 G4double bg2 = eKineticEnergy * (eTotalEnerg 285 (electron_mass_c2 * electron_ 286 287 static constexpr G4double epsfactor = 288 2. * CLHEP::electron_mass_c2 * CLHEP::elec 289 CLHEP::Bohr_radius * CLHEP::Bohr_radius / 290 G4double eps = epsfactor * bg2 / Z23; 291 292 if(eps < epsmin) 293 sigma = 2. * eps * eps; 294 else if(eps < epsmax) 295 sigma = G4Log(1. + 2. * eps) - 2. * eps / 296 else 297 sigma = G4Log(2. * eps) - 1. + 1. / eps; 298 299 sigma *= ChargeSquare * AtomicNumber * Atomi 300 301 // interpolate in AtomicNumber and beta2 302 G4double c1, c2, cc1, cc2, corr; 303 304 // get bin number in Z 305 G4int iZ = 14; 306 while((iZ >= 0) && (Zdat[iZ] >= AtomicNumber 307 iZ -= 1; 308 if(iZ == 14) 309 iZ = 13; 310 if(iZ == -1) 311 iZ = 0; 312 313 G4double ZZ1 = Zdat[iZ]; 314 G4double ZZ2 = Zdat[iZ + 1]; 315 G4double ratZ = 316 (AtomicNumber - ZZ1) * (AtomicNumber + ZZ1 317 318 static constexpr G4double Tlim = 10. * CLHEP 319 static constexpr G4double sigmafactor = 320 CLHEP::twopi * CLHEP::classic_electr_radiu 321 static const G4double beta2lim = 322 Tlim * (Tlim + 2. * CLHEP::electron_mass_c 323 ((Tlim + CLHEP::electron_mass_c2) * (Tlim 324 static const G4double bg2lim = 325 Tlim * (Tlim + 2. * CLHEP::electron_mass_c 326 (CLHEP::electron_mass_c2 * CLHEP::electron 327 328 static constexpr G4double sig0[15] = { 329 0.2672 * CLHEP::barn, 0.5922 * CLHEP::barn 330 6.235 * CLHEP::barn, 11.69 * CLHEP::barn, 331 16.12 * CLHEP::barn, 23.00 * CLHEP::barn, 332 39.95 * CLHEP::barn, 50.85 * CLHEP::barn, 333 91.15 * CLHEP::barn, 104.4 * CLHEP::barn, 334 }; 335 336 static constexpr G4double Tdat[22] = { 337 100 * CLHEP::eV, 200 * CLHEP::eV, 400 * 338 1 * CLHEP::keV, 2 * CLHEP::keV, 4 * CL 339 10 * CLHEP::keV, 20 * CLHEP::keV, 40 * C 340 100 * CLHEP::keV, 200 * CLHEP::keV, 400 * 341 1 * CLHEP::MeV, 2 * CLHEP::MeV, 4 * CL 342 10 * CLHEP::MeV, 20 * CLHEP::MeV 343 }; 344 345 if(eKineticEnergy <= Tlim) 346 { 347 // get bin number in T (beta2) 348 G4int iT = 21; 349 while((iT >= 0) && (Tdat[iT] >= eKineticEn 350 iT -= 1; 351 if(iT == 21) 352 iT = 20; 353 if(iT == -1) 354 iT = 0; 355 356 // calculate betasquare values 357 G4double T = Tdat[iT], E = T + electron_ma 358 G4double b2small = T * (E + electron_mass_ 359 360 T = Tdat[iT + 1]; 361 E = T + electron_mass_c2; 362 G4double b2big = T * (E + electron_mass_c2 363 G4double ratb2 = (beta2 - b2small) / (b2bi 364 365 if(charge < 0.) 366 { 367 c1 = celectron[iZ][iT]; 368 c2 = celectron[iZ + 1][iT]; 369 cc1 = c1 + ratZ * (c2 - c1); 370 371 c1 = celectron[iZ][iT + 1]; 372 c2 = celectron[iZ + 1][iT + 1]; 373 cc2 = c1 + ratZ * (c2 - c1); 374 375 corr = cc1 + ratb2 * (cc2 - cc1); 376 377 sigma *= sigmafactor / corr; 378 } 379 else 380 { 381 c1 = cpositron[iZ][iT]; 382 c2 = cpositron[iZ + 1][iT]; 383 cc1 = c1 + ratZ * (c2 - c1); 384 385 c1 = cpositron[iZ][iT + 1]; 386 c2 = cpositron[iZ + 1][iT + 1]; 387 cc2 = c1 + ratZ * (c2 - c1); 388 389 corr = cc1 + ratb2 * (cc2 - cc1); 390 391 sigma *= sigmafactor / corr; 392 } 393 } 394 else 395 { 396 c1 = bg2lim * sig0[iZ] * (1. + hecorr[iZ] 397 c2 = 398 bg2lim * sig0[iZ + 1] * (1. + hecorr[iZ 399 if((AtomicNumber >= ZZ1) && (AtomicNumber 400 sigma = c1 + ratZ * (c2 - c1); 401 else if(AtomicNumber < ZZ1) 402 sigma = AtomicNumber * AtomicNumber * c1 403 else if(AtomicNumber > ZZ2) 404 sigma = AtomicNumber * AtomicNumber * c2 405 } 406 return sigma; 407 } 408 409 //....oooOO0OOooo........oooOO0OOooo........oo 410 void G4UrbanAdjointMscModel::StartTracking(G4T 411 { 412 SetParticle(track->GetDynamicParticle()->Get 413 firstStep = true; 414 insideskin = false; 415 fr = facrange; 416 tlimit = tgeom = rangeinit = rangecut = geom 417 smallstep = 1.e1 418 stepmin = tlim 419 tlimitmin = 10. 420 rndmEngineMod = G4Ra 421 } 422 423 //....oooOO0OOooo........oooOO0OOooo........oo 424 G4double G4UrbanAdjointMscModel::ComputeTruePa 425 const G4Track& track, G4double& currentMinim 426 { 427 tPathLength = currentMinimal 428 const G4DynamicParticle* dp = track.GetDynam 429 430 G4StepPoint* sp = track.GetStep()->G 431 G4StepStatus stepStatus = sp->GetStepStatus( 432 couple = track.GetMaterialC 433 SetCurrentCouple(couple); 434 currentMaterialIndex = couple->GetIndex(); 435 currentKinEnergy = dp->GetKineticEnergy( 436 437 currentRange = GetRange(particle, currentKin 438 lambda0 = GetTransportMeanFreePath(part 439 tPathLength = min(tPathLength, currentRange 440 441 // set flag to default values 442 Zeff = couple->GetMaterial()->GetIonisation( 443 444 if(Zold != Zeff) 445 UpdateCache(); 446 447 // stop here if small step 448 if(tPathLength < tlimitminfix) 449 { 450 latDisplasment = false; 451 return ConvertTrueToGeom(tPathLength, curr 452 } 453 454 // upper limit for the straight line distanc 455 // for electrons and positrons 456 G4double distance = currentRange; 457 // for muons, hadrons 458 if(mass > masslimite) 459 { 460 distance *= (1.15 - 9.76e-4 * Zeff); 461 } 462 else 463 { 464 distance *= (1.20 - Zeff * (1.62e-2 - 9.22 465 } 466 presafety = sp->GetSafety(); 467 468 // far from geometry boundary 469 if(distance < presafety) 470 { 471 latDisplasment = false; 472 return ConvertTrueToGeom(tPathLength, curr 473 } 474 475 latDisplasment = latDispla 476 static constexpr G4double invmev = 1.0 / CLH 477 478 // standard version 479 if(steppingAlgorithm == fUseDistanceToBounda 480 { 481 // compute geomlimit and presafety 482 geomlimit = ComputeGeomLimit(track, presaf 483 484 // is it far from boundary ? 485 if(distance < presafety) 486 { 487 latDisplasment = false; 488 return ConvertTrueToGeom(tPathLength, cu 489 } 490 491 smallstep += 1.; 492 insideskin = false; 493 494 // initialisation at firs step and at the 495 if(firstStep || (stepStatus == fGeomBounda 496 { 497 rangeinit = currentRange; 498 if(!firstStep) 499 { 500 smallstep = 1.; 501 } 502 503 // define stepmin here (it depends on la 504 // rough estimation of lambda_elastic/la 505 G4double rat = currentKinEnergy * invmev 506 rat = 1.e-3 / (rat * (10. + rat 507 // stepmin ~ lambda_elastic 508 stepmin = rat * lambda0; 509 skindepth = skin * stepmin; 510 tlimitmin = max(10 * stepmin, tlimitminf 511 512 // constraint from the geometry 513 if((geomlimit < geombig) && (geomlimit > 514 { 515 // geomlimit is a geometrical step len 516 // transform it to true path length (e 517 if((1. - geomlimit / lambda0) > 0.) 518 geomlimit = -lambda0 * G4Log(1. - ge 519 520 if(stepStatus == fGeomBoundary) 521 tgeom = geomlimit / facgeom; 522 else 523 tgeom = 2. * geomlimit / facgeom; 524 } 525 else 526 tgeom = geombig; 527 } 528 529 // step limit 530 tlimit = facrange * rangeinit; 531 532 // lower limit for tlimit 533 tlimit = max(tlimit, tlimitmin); 534 tlimit = min(tlimit, tgeom); 535 536 // shortcut 537 if((tPathLength < tlimit) && (tPathLength 538 (smallstep > skin) && (tPathLength < ge 539 { 540 return ConvertTrueToGeom(tPathLength, cu 541 } 542 543 // step reduction near to boundary 544 if(smallstep <= skin) 545 { 546 tlimit = stepmin; 547 insideskin = true; 548 } 549 else if(geomlimit < geombig) 550 { 551 if(geomlimit > skindepth) 552 { 553 tlimit = min(tlimit, geomlimit - 0.999 554 } 555 else 556 { 557 insideskin = true; 558 tlimit = min(tlimit, stepmin); 559 } 560 } 561 562 tlimit = max(tlimit, stepmin); 563 564 // randomise if not 'small' step and step 565 if((tlimit < tPathLength) && (smallstep > 566 { 567 tPathLength = min(tPathLength, Randomize 568 } 569 else 570 { 571 tPathLength = min(tPathLength, tlimit); 572 } 573 } 574 // for 'normal' simulation with or without m 575 // there no small step/single scattering at 576 else if(steppingAlgorithm == fUseSafety) 577 { 578 if(stepStatus != fGeomBoundary) 579 { 580 presafety = ComputeSafety(sp->GetPositio 581 } 582 // is far from boundary 583 if(distance < presafety) 584 { 585 latDisplasment = false; 586 return ConvertTrueToGeom(tPathLength, cu 587 } 588 589 if(firstStep || (stepStatus == fGeomBounda 590 { 591 rangeinit = currentRange; 592 fr = facrange; 593 // Geant4 version 9.1 like stepping for 594 if(mass < masslimite) 595 { 596 rangeinit = max(rangeinit, lambda0); 597 if(lambda0 > lambdalimit) 598 { 599 fr *= (0.75 + 0.25 * lambda0 / lambd 600 } 601 } 602 // lower limit for tlimit 603 G4double rat = currentKinEnergy * invmev 604 rat = 1.e-3 / (rat * (10 + rat) 605 stepmin = lambda0 * rat; 606 tlimitmin = max(10 * stepmin, tlimitm 607 } 608 609 // step limit 610 tlimit = max(fr * rangeinit, facsafety * p 611 612 // lower limit for tlimit 613 tlimit = max(tlimit, tlimitmin); 614 615 // randomise if step determined by msc 616 if(tlimit < tPathLength) 617 { 618 tPathLength = min(tPathLength, Randomize 619 } 620 else 621 { 622 tPathLength = min(tPathLength, tlimit); 623 } 624 } 625 // new stepping mode UseSafetyPlus 626 else if(steppingAlgorithm == fUseSafetyPlus) 627 { 628 if(stepStatus != fGeomBoundary) 629 { 630 presafety = ComputeSafety(sp->GetPositio 631 } 632 633 // is far from boundary 634 if(distance < presafety) 635 { 636 latDisplasment = false; 637 return ConvertTrueToGeom(tPathLength, cu 638 } 639 640 if(firstStep || (stepStatus == fGeomBounda 641 { 642 rangeinit = currentRange; 643 fr = facrange; 644 rangecut = geombig; 645 if(mass < masslimite) 646 { 647 G4int index = 1; 648 if(charge > 0.) 649 index = 2; 650 rangecut = couple->GetProductionCuts() 651 if(lambda0 > lambdalimit) 652 { 653 fr *= (0.84 + 0.16 * lambda0 / lambd 654 } 655 } 656 // lower limit for tlimit 657 G4double rat = currentKinEnergy * invmev 658 rat = 1.e-3 / (rat * (10 + rat) 659 stepmin = lambda0 * rat; 660 tlimitmin = max(10 * stepmin, tlimitm 661 } 662 // step limit 663 tlimit = max(fr * rangeinit, facsafety * p 664 665 // lower limit for tlimit 666 tlimit = max(tlimit, tlimitmin); 667 668 // condition for tPathLength from drr and 669 if(currentRange > finalr) 670 { 671 G4double tmax = 672 drr * currentRange + finalr * (1. - dr 673 tPathLength = min(tPathLength, tmax); 674 } 675 676 // condition safety 677 if(currentRange > rangecut) 678 { 679 if(firstStep) 680 { 681 tPathLength = min(tPathLength, facsafe 682 } 683 else if(stepStatus != fGeomBoundary && p 684 { 685 tPathLength = min(tPathLength, presafe 686 } 687 } 688 689 // randomise if step determined by msc 690 if(tPathLength < tlimit) 691 { 692 tPathLength = min(tPathLength, Randomize 693 } 694 else 695 { 696 tPathLength = min(tPathLength, tlimit); 697 } 698 } 699 700 // version similar to 7.1 (needed for some e 701 else 702 { 703 if(stepStatus == fGeomBoundary) 704 { 705 if(currentRange > lambda0) 706 { 707 tlimit = facrange * currentRange; 708 } 709 else 710 { 711 tlimit = facrange * lambda0; 712 } 713 714 tlimit = max(tlimit, tlimitmin); 715 } 716 // randomise if step determined by msc 717 if(tlimit < tPathLength) 718 { 719 tPathLength = min(tPathLength, Randomize 720 } 721 else 722 { 723 tPathLength = min(tPathLength, tlimit); 724 } 725 } 726 firstStep = false; 727 return ConvertTrueToGeom(tPathLength, curren 728 } 729 730 //....oooOO0OOooo........oooOO0OOooo........oo 731 G4double G4UrbanAdjointMscModel::ComputeGeomPa 732 { 733 lambdaeff = lambda0; 734 par1 = -1.; 735 par2 = par3 = 0.; 736 737 // this correction needed to run MSC with eI 738 // and makes no harm for a normal run 739 tPathLength = std::min(tPathLength, currentR 740 741 // do the true -> geom transformation 742 zPathLength = tPathLength; 743 744 // z = t for very small tPathLength 745 if(tPathLength < tlimitminfix2) 746 return zPathLength; 747 748 G4double tau = tPathLength / lambda0; 749 750 if((tau <= tausmall) || insideskin) 751 { 752 zPathLength = min(tPathLength, lambda0); 753 } 754 else if(tPathLength < currentRange * dtrl) 755 { 756 if(tau < taulim) 757 zPathLength = tPathLength * (1. - 0.5 * 758 else 759 zPathLength = lambda0 * (1. - G4Exp(-tau 760 } 761 else if(currentKinEnergy < mass || tPathLeng 762 { 763 par1 = 1. / currentRange; 764 par2 = 1. / (par1 * lambda0); 765 par3 = 1. + par2; 766 if(tPathLength < currentRange) 767 { 768 zPathLength = 769 (1. - G4Exp(par3 * G4Log(1. - tPathLen 770 (par1 * par3); 771 } 772 else 773 { 774 zPathLength = 1. / (par1 * par3); 775 } 776 } 777 else 778 { 779 G4double rfin = max(currentRange - tPat 780 G4double T1 = GetEnergy(particle, rfi 781 G4double lambda1 = GetTransportMeanFreePat 782 783 par1 = (lambda0 - lambda1) / (lambd 784 par2 = 1. / (par1 * lambda0); 785 par3 = 1. + par2; 786 zPathLength = (1. - G4Exp(par3 * G4Log(lam 787 } 788 789 zPathLength = min(zPathLength, lambda0); 790 return zPathLength; 791 } 792 793 //....oooOO0OOooo........oooOO0OOooo........oo 794 G4double G4UrbanAdjointMscModel::ComputeTrueSt 795 { 796 // step defined other than transportation 797 if(geomStepLength == zPathLength) 798 { 799 return tPathLength; 800 } 801 802 zPathLength = geomStepLength; 803 804 // t = z for very small step 805 if(geomStepLength < tlimitminfix2) 806 { 807 tPathLength = geomStepLength; 808 809 // recalculation 810 } 811 else 812 { 813 G4double tlength = geomStepLength; 814 if((geomStepLength > lambda0 * tausmall) & 815 { 816 if(par1 < 0.) 817 { 818 tlength = -lambda0 * G4Log(1. - geomSt 819 } 820 else 821 { 822 if(par1 * par3 * geomStepLength < 1.) 823 { 824 tlength = 825 (1. - G4Exp(G4Log(1. - par1 * par3 826 par1; 827 } 828 else 829 { 830 tlength = currentRange; 831 } 832 } 833 834 if(tlength < geomStepLength) 835 { 836 tlength = geomStepLength; 837 } 838 else if(tlength > tPathLength) 839 { 840 tlength = tPathLength; 841 } 842 } 843 tPathLength = tlength; 844 } 845 846 return tPathLength; 847 } 848 849 //....oooOO0OOooo........oooOO0OOooo........oo 850 G4ThreeVector& G4UrbanAdjointMscModel::SampleS 851 const G4ThreeVector& oldDirection, G4double 852 { 853 fDisplacement.set(0.0, 0.0, 0.0); 854 G4double kineticEnergy = currentKinEnergy; 855 if(tPathLength > currentRange * dtrl) 856 { 857 kineticEnergy = GetEnergy(particle, curren 858 } 859 else 860 { 861 kineticEnergy -= tPathLength * GetDEDX(par 862 } 863 864 if((kineticEnergy <= eV) || (tPathLength <= 865 (tPathLength < tausmall * lambda0)) 866 { 867 return fDisplacement; 868 } 869 870 G4double cth = SampleCosineTheta(tPathLength 871 872 // protection against 'bad' cth values 873 if(std::fabs(cth) >= 1.0) 874 { 875 return fDisplacement; 876 } 877 878 G4double sth = sqrt((1.0 - cth) * (1.0 + ct 879 G4double phi = twopi * rndmEngineMod->flat( 880 G4double dirx = sth * cos(phi); 881 G4double diry = sth * sin(phi); 882 883 G4ThreeVector newDirection(dirx, diry, cth); 884 newDirection.rotateUz(oldDirection); 885 886 fParticleChange->ProposeMomentumDirection(ne 887 888 if(latDisplasment && currentTau >= tausmall) 889 { 890 if(displacementFlag) 891 { 892 SampleDisplacementNew(cth, phi); 893 } 894 else 895 { 896 SampleDisplacement(sth, phi); 897 } 898 fDisplacement.rotateUz(oldDirection); 899 } 900 return fDisplacement; 901 } 902 903 //....oooOO0OOooo........oooOO0OOooo........oo 904 G4double G4UrbanAdjointMscModel::SampleCosineT 905 906 { 907 G4double cth = 1.; 908 G4double tau = trueStepLength / lambda0; 909 currentTau = tau; 910 lambdaeff = lambda0; 911 912 G4double lambda1 = GetTransportMeanFreePath( 913 if(std::fabs(lambda1 - lambda0) > lambda0 * 914 { 915 // mean tau value 916 tau = trueStepLength * G4Log(lambda0 / lam 917 } 918 919 currentTau = tau; 920 lambdaeff = trueStepLength / currentT 921 currentRadLength = couple->GetMaterial()->Ge 922 923 if(tau >= taubig) 924 { 925 cth = -1. + 2. * rndmEngineMod->flat(); 926 } 927 else if(tau >= tausmall) 928 { 929 static const G4double numlim = 0.01; 930 G4double xmeanth, x2meanth; 931 if(tau < numlim) 932 { 933 xmeanth = 1.0 - tau * (1.0 - 0.5 * tau) 934 x2meanth = 1.0 - tau * (5.0 - 6.25 * tau 935 } 936 else 937 { 938 xmeanth = G4Exp(-tau); 939 x2meanth = (1. + 2. * G4Exp(-2.5 * tau)) 940 } 941 942 // too large step of low-energy particle 943 G4double relloss = 1. - Ki 944 static const G4double rellossmax = 0.50; 945 if(relloss > rellossmax) 946 { 947 return SimpleScattering(xmeanth, x2meant 948 } 949 // is step extreme small ? 950 G4bool extremesmallstep = false; 951 G4double tsmall = std::min(tlimitm 952 G4double theta0 = 0.; 953 if(trueStepLength > tsmall) 954 { 955 theta0 = ComputeTheta0(trueStepLength, K 956 } 957 else 958 { 959 theta0 = 960 sqrt(trueStepLength / tsmall) * Comput 961 extremesmallstep = true; 962 } 963 964 static constexpr G4double theta0max = CLHE 965 966 // protection for very small angles 967 G4double theta2 = theta0 * theta0; 968 969 if(theta2 < tausmall) 970 { 971 return cth; 972 } 973 974 if(theta0 > theta0max) 975 { 976 return SimpleScattering(xmeanth, x2meant 977 } 978 979 G4double x = theta2 * (1.0 - theta2 / 12.) 980 if(theta2 > numlim) 981 { 982 G4double sth = 2 * sin(0.5 * theta0); 983 x = sth * sth; 984 } 985 986 // parameter for tail 987 G4double ltau = G4Log(tau); 988 G4double u = G4Exp(ltau / 6.); 989 if(extremesmallstep) 990 u = G4Exp(G4Log(tsmall / lambda0) / 6.); 991 G4double xx = G4Log(lambdaeff / currentRa 992 G4double xsi = coeffc1 + u * (coeffc2 + co 993 994 // tail should not be too big 995 if(xsi < 1.9) 996 { 997 xsi = 1.9; 998 } 999 1000 G4double c = xsi; 1001 1002 if(std::abs(c - 3.) < 0.001) 1003 { 1004 c = 3.001; 1005 } 1006 else if(std::abs(c - 2.) < 0.001) 1007 { 1008 c = 2.001; 1009 } 1010 1011 G4double c1 = c - 1.; 1012 1013 G4double ea = G4Exp(-xsi); 1014 G4double eaa = 1. - ea; 1015 G4double xmean1 = 1. - (1. - (1. + xsi) * 1016 G4double x0 = 1. - xsi * x; 1017 1018 if(xmean1 <= 0.999 * xmeanth) 1019 { 1020 return SimpleScattering(xmeanth, x2mean 1021 } 1022 // from continuity of derivatives 1023 G4double b = 1. + (c - xsi) * x; 1024 1025 G4double b1 = b + 1.; 1026 G4double bx = c * x; 1027 1028 G4double eb1 = G4Exp(G4Log(b1) * c1); 1029 G4double ebx = G4Exp(G4Log(bx) * c1); 1030 G4double d = ebx / eb1; 1031 1032 G4double xmean2 = (x0 + d - (bx - b1 * d) 1033 1034 G4double f1x0 = ea / eaa; 1035 G4double f2x0 = c1 / (c * (1. - d)); 1036 G4double prob = f2x0 / (f1x0 + f2x0); 1037 1038 G4double qprob = xmeanth / (prob * xmean1 1039 1040 if(rndmEngineMod->flat() < qprob) 1041 { 1042 G4double var = 0; 1043 if(rndmEngineMod->flat() < prob) 1044 { 1045 cth = 1. + G4Log(ea + rndmEngineMod-> 1046 } 1047 else 1048 { 1049 var = (1.0 - d) * rndmEngineMod->flat 1050 if(var < numlim * d) 1051 { 1052 var /= (d * c1); 1053 cth = -1.0 + var * (1.0 - 0.5 * var 1054 } 1055 else 1056 { 1057 cth = 1. + x * (c - xsi - c * G4Exp 1058 } 1059 } 1060 } 1061 else 1062 { 1063 cth = -1. + 2. * rndmEngineMod->flat(); 1064 } 1065 } 1066 return cth; 1067 } 1068 1069 //....oooOO0OOooo........oooOO0OOooo........o 1070 G4double G4UrbanAdjointMscModel::ComputeTheta 1071 1072 { 1073 // for all particles take the width of the 1074 // from a parametrization similar to the H 1075 // ( Highland formula: Particle Physics Boo 1076 G4double invbetacp = 1077 std::sqrt((currentKinEnergy + mass) * (Ki 1078 (currentKinEnergy * (currentKin 1079 KineticEnergy * (KineticEnergy 1080 G4double y = trueStepLength / currentRadLen 1081 1082 if(particle == positron) 1083 { 1084 static constexpr G4double xl = 0.6; 1085 static constexpr G4double xh = 0.9; 1086 static constexpr G4double e = 113.0; 1087 G4double corr; 1088 1089 G4double tau = std::sqrt(currentKinEnergy 1090 G4double x = std::sqrt(tau * (tau + 2.) 1091 G4double a = 0.994 - 4.08e-3 * Zeff; 1092 G4double b = 7.16 + (52.6 + 365. / Zeff 1093 G4double c = 1.000 - 4.47e-3 * Zeff; 1094 G4double d = 1.21e-3 * Zeff; 1095 if(x < xl) 1096 { 1097 corr = a * (1. - G4Exp(-b * x)); 1098 } 1099 else if(x > xh) 1100 { 1101 corr = c + d * G4Exp(e * (x - 1.)); 1102 } 1103 else 1104 { 1105 G4double yl = a * (1. - G4Exp(-b * xl)) 1106 G4double yh = c + d * G4Exp(e * (xh - 1 1107 G4double y0 = (yh - yl) / (xh - xl); 1108 G4double y1 = yl - y0 * xl; 1109 corr = y0 * x + y1; 1110 } 1111 y *= corr * (1. + Zeff * (1.84035e-4 * Ze 1112 } 1113 1114 static constexpr G4double c_highland = 13.6 1115 G4double theta0 = c_highland * std::abs(cha 1116 1117 // correction factor from e- scattering dat 1118 theta0 *= (coeffth1 + coeffth2 * G4Log(y)); 1119 return theta0; 1120 } 1121 1122 //....oooOO0OOooo........oooOO0OOooo........o 1123 void G4UrbanAdjointMscModel::SampleDisplaceme 1124 { 1125 G4double rmax = 1126 sqrt((tPathLength - zPathLength) * (tPath 1127 1128 static constexpr G4double third = 1. / 3.; 1129 G4double r = rmax * G4Exp(G4Log(rndmEngineM 1130 1131 if(r > 0.) 1132 { 1133 static constexpr G4double kappa = 2.5; 1134 static constexpr G4double kappami1 = 1.5; 1135 1136 G4double latcorr = 0.; 1137 if((currentTau >= tausmall) && !insideski 1138 { 1139 if(currentTau < taulim) 1140 { 1141 latcorr = lambdaeff * kappa * current 1142 (1. - (kappa + 1.) * curren 1143 } 1144 else 1145 { 1146 G4double etau = 0.; 1147 if(currentTau < taubig) 1148 { 1149 etau = G4Exp(-currentTau); 1150 } 1151 latcorr = -kappa * currentTau; 1152 latcorr = G4Exp(latcorr) / kappami1; 1153 latcorr += 1. - kappa * etau / kappam 1154 latcorr *= 2. * lambdaeff * third; 1155 } 1156 } 1157 latcorr = std::min(latcorr, r); 1158 1159 // sample direction of lateral displaceme 1160 // compute it from the lateral correlatio 1161 G4double Phi = 0.; 1162 if(std::abs(r * sth) < latcorr) 1163 { 1164 Phi = twopi * rndmEngineMod->flat(); 1165 } 1166 else 1167 { 1168 G4double psi = std::acos(latcorr / (r * 1169 if(rndmEngineMod->flat() < 0.5) 1170 { 1171 Phi = phi + psi; 1172 } 1173 else 1174 { 1175 Phi = phi - psi; 1176 } 1177 } 1178 fDisplacement.set(r * std::cos(Phi), r * 1179 } 1180 } 1181 1182 //....oooOO0OOooo........oooOO0OOooo........o 1183 void G4UrbanAdjointMscModel::SampleDisplaceme 1184 { 1185 // sample displacement r 1186 1187 G4double rmax = 1188 sqrt((tPathLength - zPathLength) * (tPath 1189 // u = (r/rmax)**2 , v=1-u 1190 // paramerization from ss simulation 1191 // f(u) = p0*exp(p1*log(v)-p2*v)+v*(p3+p4*v 1192 G4double u, v, rej; 1193 G4int count = 0; 1194 1195 static constexpr G4double reps = 1.e-6; 1196 static constexpr G4double rp0 = 2.2747e+4; 1197 static constexpr G4double rp1 = 4.5980e+0; 1198 static constexpr G4double rp2 = 1.5580e+1; 1199 static constexpr G4double rp3 = 7.1287e-1; 1200 static constexpr G4double rp4 = -5.7069e-1 1201 1202 do 1203 { 1204 u = reps + (1. - 2. * reps) * rndmEngin 1205 v = 1. - u; 1206 rej = rp0 * G4Exp(rp1 * G4Log(v) - rp2 * 1207 } while(rndmEngineMod->flat() > rej && ++co 1208 G4double r = rmax * sqrt(u); 1209 1210 if(r > 0.) 1211 { 1212 // sample Phi using lateral correlation 1213 // v = Phi-phi = acos(latcorr/(r*sth)) 1214 // v has a universal distribution which c 1215 // simulation as 1216 // f(v) = 1.49e-2*exp(-v**2/(2*0.320))+2. 1217 // 1.96e-5*exp(8.42e-1*log(1.+1.45 1218 static constexpr G4double probv1 = 0.3055 1219 static constexpr G4double probv2 = 0.9551 1220 static constexpr G4double vhigh = 3.15; 1221 static const G4double w2v = 1. / G4Exp(30 1222 static const G4double w3v = 1. / G4Exp(-1 1223 1224 G4double Phi; 1225 G4double random = rndmEngineMod->flat(); 1226 if(random < probv1) 1227 { 1228 do 1229 { 1230 v = G4RandGauss::shoot(rndmEngineMod, 1231 } while(std::abs(v) >= vhigh); 1232 Phi = phi + v; 1233 } 1234 else 1235 { 1236 if(random < probv2) 1237 { 1238 v = (-1. + 1239 1. / G4Exp(G4Log(1. - rndmEngine 1240 6.30e-2; 1241 } 1242 else 1243 { 1244 v = (-1. + 1. / G4Exp(G4Log(1. - rndm 1245 -1.842)) / 1246 1.45e1; 1247 } 1248 1249 random = rndmEngineMod->flat(); 1250 if(random < 0.5) 1251 { 1252 Phi = phi + v; 1253 } 1254 else 1255 { 1256 Phi = phi - v; 1257 } 1258 } 1259 fDisplacement.set(r * std::cos(Phi), r * 1260 } 1261 } 1262 1263 //....oooOO0OOooo........oooOO0OOooo........o 1264