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 // GEANT4 Class file 29 // 30 // 31 // File name: G4UrbanMscModel 32 // 33 // Author: Laszlo Urban 34 // 35 // Creation date: 19.02.2013 36 // 37 // Created from G4UrbanMscModel96 38 // 39 // New parametrization for theta0 40 // Correction for very small step length 41 // 42 // Class Description: 43 // 44 // Implementation of the model of multiple sca 45 // H.W.Lewis Phys Rev 78 (1950) 526 and others 46 47 // ------------------------------------------- 48 // In its present form the model can be used 49 // of the e-/e+ multiple scattering 50 51 52 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oo 54 55 #include "G4UrbanMscModel.hh" 56 #include "G4PhysicalConstants.hh" 57 #include "G4SystemOfUnits.hh" 58 #include "Randomize.hh" 59 #include "G4Positron.hh" 60 #include "G4EmParameters.hh" 61 #include "G4ParticleChangeForMSC.hh" 62 #include "G4ProductionCutsTable.hh" 63 64 #include "G4Poisson.hh" 65 #include "G4Pow.hh" 66 #include "G4Log.hh" 67 #include "G4Exp.hh" 68 #include "G4AutoLock.hh" 69 70 //....oooOO0OOooo........oooOO0OOooo........oo 71 72 std::vector<G4UrbanMscModel::mscData*> G4Urban 73 74 namespace 75 { 76 G4Mutex theUrbanMutex = G4MUTEX_INITIALIZER; 77 } 78 79 //....oooOO0OOooo........oooOO0OOooo........oo 80 81 G4UrbanMscModel::G4UrbanMscModel(const G4Strin 82 : G4VMscModel(nam) 83 { 84 masslimite = 0.6*CLHEP::MeV; 85 fr = 0.02; 86 taubig = 8.0; 87 tausmall = 1.e-16; 88 taulim = 1.e-6; 89 currentTau = taulim; 90 tlimitminfix = 0.01*CLHEP::nm; 91 tlimitminfix2 = 1.*CLHEP::nm; 92 stepmin = tlimitminfix; 93 smallstep = 1.e10; 94 currentRange = 0.; 95 rangeinit = 0.; 96 tlimit = 1.e10*CLHEP::mm; 97 tlimitmin = 10.*tlimitminfix; 98 tgeom = 1.e50*CLHEP::mm; 99 geombig = tgeom; 100 geommin = 1.e-3*CLHEP::mm; 101 geomlimit = geombig; 102 presafety = 0.; 103 104 positron = G4Positron::Positron(); 105 rndmEngineMod = G4Random::getTheEngine(); 106 107 drr = 0.35; 108 finalr = 10.*CLHEP::um; 109 110 tlow = 5.*CLHEP::keV; 111 invmev = 1.0/CLHEP::MeV; 112 113 skindepth = skin*stepmin; 114 115 mass = CLHEP::proton_mass_c2; 116 charge = chargeSquare = 1.0; 117 currentKinEnergy = currentRadLength = lambda 118 = zPathLength = par1 = par2 = par3 = rndma 119 currentLogKinEnergy = LOG_EKIN_MIN; 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oo 123 124 G4UrbanMscModel::~G4UrbanMscModel() 125 { 126 if(isFirstInstance) { 127 for(auto const & ptr : msc) { delete ptr; 128 msc.clear(); 129 } 130 } 131 132 //....oooOO0OOooo........oooOO0OOooo........oo 133 134 void G4UrbanMscModel::Initialise(const G4Parti 135 const G4DataV 136 { 137 // set values of some data members 138 SetParticle(p); 139 fParticleChange = GetParticleChangeForMSC(p) 140 InitialiseParameters(p); 141 142 latDisplasmentbackup = latDisplasment; 143 144 // if model is locked parameters should be d 145 if(!IsLocked()) { 146 dispAlg96 = G4EmParameters::Instance()->La 147 fPosiCorrection = G4EmParameters::Instance 148 } 149 150 // initialise cache only once 151 if(0 == msc.size()) { 152 G4AutoLock l(&theUrbanMutex); 153 if(0 == msc.size()) { 154 isFirstInstance = true; 155 msc.resize(1, nullptr); 156 } 157 l.unlock(); 158 } 159 // initialise cache for each new run 160 if(isFirstInstance) { InitialiseModelCache() 161 162 /* 163 G4cout << "### G4UrbanMscModel::Initialise d 164 << p->GetParticleName() << " type= " << ste 165 G4cout << " RangeFact= " << facrange << " 166 << " SafetyFact= " << facsafety << " Lambda 167 << G4endl; 168 */ 169 } 170 171 //....oooOO0OOooo........oooOO0OOooo........oo 172 173 G4double G4UrbanMscModel::ComputeCrossSectionP 174 const G4ParticleD 175 G4double ki 176 G4double at 177 G4double, G 178 { 179 static const G4double epsmin = 1.e-4 , epsma 180 181 static const G4double Zdat[15] = { 4., 6., 182 50., 56., 183 184 // corr. factors for e-/e+ lambda for T <= T 185 static const G4double celectron[15][22] = 186 {{1.125,1.072,1.051,1.047,1.047,1.05 187 1.054,1.057,1.062,1.069,1.075,1.09 188 1.112,1.108,1.100,1.093,1.089,1.08 189 {1.408,1.246,1.143,1.096,1.077,1.05 190 1.052,1.053,1.058,1.065,1.072,1.08 191 1.109,1.105,1.097,1.090,1.086,1.08 192 {2.833,2.268,1.861,1.612,1.486,1.30 193 1.136,1.114,1.106,1.106,1.109,1.11 194 1.131,1.124,1.113,1.104,1.099,1.09 195 {3.879,3.016,2.380,2.007,1.818,1.53 196 1.190,1.133,1.107,1.099,1.098,1.10 197 1.112,1.105,1.096,1.089,1.085,1.09 198 {6.937,4.330,2.886,2.256,1.987,1.62 199 1.203,1.122,1.080,1.065,1.061,1.06 200 1.073,1.070,1.064,1.059,1.056,1.05 201 {9.616,5.708,3.424,2.551,2.204,1.76 202 1.256,1.155,1.099,1.077,1.070,1.06 203 1.074,1.070,1.063,1.059,1.056,1.05 204 {11.72,6.364,3.811,2.806,2.401,1.88 205 1.300,1.180,1.112,1.082,1.073,1.06 206 1.068,1.064,1.059,1.054,1.051,1.05 207 {18.08,8.601,4.569,3.183,2.662,2.02 208 1.339,1.195,1.108,1.068,1.053,1.04 209 1.039,1.037,1.034,1.031,1.030,1.03 210 {18.22,10.48,5.333,3.713,3.115,2.36 211 1.498,1.301,1.171,1.105,1.077,1.04 212 1.031,1.028,1.024,1.022,1.021,1.02 213 {14.14,10.65,5.710,3.929,3.266,2.45 214 1.528,1.319,1.178,1.106,1.075,1.04 215 1.020,1.017,1.015,1.013,1.013,1.02 216 {14.11,11.73,6.312,4.240,3.478,2.56 217 1.569,1.342,1.186,1.102,1.065,1.02 218 0.995,0.993,0.993,0.993,0.993,1.01 219 {22.76,20.01,8.835,5.287,4.144,2.90 220 1.677,1.410,1.224,1.121,1.073,1.01 221 0.974,0.972,0.973,0.974,0.975,0.98 222 {50.77,40.85,14.13,7.184,5.284,3.43 223 1.837,1.512,1.283,1.153,1.091,1.01 224 0.950,0.947,0.949,0.952,0.954,0.96 225 {65.87,59.06,15.87,7.570,5.567,3.65 226 1.939,1.579,1.325,1.178,1.108,1.01 227 0.941,0.938,0.940,0.944,0.946,0.95 228 {55.60,47.34,15.92,7.810,5.755,3.76 229 1.985,1.609,1.343,1.188,1.113,1.01 230 0.933,0.930,0.933,0.936,0.939,0.94 231 232 static const G4double cpositron[15][22] = { 233 {2.589,2.044,1.658,1.446,1.347,1.21 234 1.097,1.083,1.080,1.086,1.092,1.10 235 1.131,1.126,1.117,1.108,1.103,1.10 236 {3.904,2.794,2.079,1.710,1.543,1.32 237 1.122,1.096,1.089,1.092,1.098,1.11 238 1.138,1.132,1.122,1.113,1.108,1.10 239 {7.970,6.080,4.442,3.398,2.872,2.12 240 1.357,1.246,1.194,1.179,1.178,1.18 241 1.203,1.190,1.173,1.159,1.151,1.14 242 {9.714,7.607,5.747,4.493,3.815,2.77 243 1.553,1.353,1.253,1.219,1.211,1.21 244 1.225,1.210,1.191,1.175,1.166,1.17 245 {17.97,12.95,8.628,6.065,4.849,3.22 246 1.624,1.382,1.259,1.214,1.202,1.20 247 1.217,1.203,1.184,1.169,1.160,1.15 248 {24.83,17.06,10.84,7.355,5.767,3.70 249 1.759,1.465,1.311,1.252,1.234,1.22 250 1.237,1.222,1.201,1.184,1.174,1.15 251 {23.26,17.15,11.52,8.049,6.375,4.11 252 1.880,1.535,1.353,1.281,1.258,1.24 253 1.252,1.234,1.212,1.194,1.183,1.17 254 {22.33,18.01,12.86,9.212,7.336,4.70 255 2.015,1.602,1.385,1.297,1.268,1.25 256 1.254,1.237,1.214,1.195,1.185,1.17 257 {33.91,24.13,15.71,10.80,8.507,5.46 258 2.407,1.873,1.564,1.425,1.374,1.33 259 1.312,1.288,1.258,1.235,1.221,1.20 260 {32.14,24.11,16.30,11.40,9.015,5.78 261 2.490,1.925,1.596,1.447,1.391,1.34 262 1.320,1.294,1.264,1.240,1.226,1.21 263 {29.51,24.07,17.19,12.28,9.766,6.23 264 2.602,1.995,1.641,1.477,1.414,1.35 265 1.328,1.302,1.270,1.245,1.231,1.23 266 {38.19,30.85,21.76,15.35,12.07,7.52 267 2.926,2.188,1.763,1.563,1.484,1.40 268 1.361,1.330,1.294,1.267,1.251,1.23 269 {49.71,39.80,27.96,19.63,15.36,9.40 270 3.417,2.478,1.944,1.692,1.589,1.48 271 1.409,1.372,1.330,1.298,1.280,1.25 272 {59.25,45.08,30.36,20.83,16.15,9.83 273 3.641,2.648,2.064,1.779,1.661,1.53 274 1.442,1.400,1.354,1.319,1.299,1.27 275 {56.38,44.29,30.50,21.18,16.51,10.1 276 3.752,2.724,2.116,1.817,1.692,1.55 277 1.456,1.412,1.364,1.328,1.307,1.28 278 279 //data/corrections for T > Tlim 280 281 static const G4double hecorr[15] = { 282 120.70, 117.50, 105.00, 92.92, 79.23, 74. 283 57.39, 41.97, 36.14, 24.53, 10.21, -7.8 284 -22.30}; 285 286 G4double sigma; 287 SetParticle(part); 288 289 G4double Z23 = G4Pow::GetInstance()->Z23(G4l 290 291 // correction if particle .ne. e-/e+ 292 // compute equivalent kinetic energy 293 // lambda depends on p*beta .... 294 295 G4double eKineticEnergy = kinEnergy; 296 297 if(mass > CLHEP::electron_mass_c2) 298 { 299 G4double TAU = kinEnergy/mass ; 300 G4double c = mass*TAU*(TAU+2.)/(CLHEP::el 301 G4double w = c-2.; 302 G4double tau = 0.5*(w+std::sqrt(w*w+4.*c) 303 eKineticEnergy = CLHEP::electron_mass_c2* 304 } 305 306 G4double eTotalEnergy = eKineticEnergy + CLH 307 G4double beta2 = eKineticEnergy*(eTotalEnerg 308 /(eTotalEnerg 309 G4double bg2 = eKineticEnergy*(eTotalEnerg 310 /(CLHEP::elec 311 312 static const G4double epsfactor = 2.*CLHEP:: 313 CLHEP::electron_mass_c2*CLHEP::Bohr_radius 314 /(CLHEP::hbarc*CLHEP::hbarc); 315 G4double eps = epsfactor*bg2/Z23; 316 317 if (eps<epsmin) sigma = 2.*eps*eps; 318 else if(eps<epsmax) sigma = G4Log(1.+2.*eps 319 else sigma = G4Log(2.*eps)-1 320 321 sigma *= chargeSquare*atomicNumber*atomicNum 322 323 // interpolate in AtomicNumber and beta2 324 G4double c1,c2,cc1; 325 326 // get bin number in Z 327 G4int iZ = 14; 328 // Loop checking, 03-Aug-2015, Vladimir Ivan 329 while ((iZ>=0)&&(Zdat[iZ]>=atomicNumber)) { 330 331 iZ = std::min(std::max(iZ, 0), 13); 332 333 G4double ZZ1 = Zdat[iZ]; 334 G4double ZZ2 = Zdat[iZ+1]; 335 G4double ratZ = (atomicNumber-ZZ1)*(atomicNu 336 ((ZZ2-ZZ1)*(ZZ2+ZZ1)); 337 338 static const G4double Tlim = 10.*CLHEP::MeV; 339 static const G4double sigmafactor = 340 CLHEP::twopi*CLHEP::classic_electr_radius* 341 static const G4double beta2lim = Tlim*(Tlim+ 342 ((Tlim+CLHEP::electron_mass_c2)*(Tlim+CLHE 343 static const G4double bg2lim = Tlim*(Tlim+ 344 (CLHEP::electron_mass_c2*CLHEP::electron_m 345 346 static const G4double sig0[15] = { 347 0.2672*CLHEP::barn, 0.5922*CLHEP::barn, 348 11.69*CLHEP::barn , 13.24*CLHEP::barn , 349 35.13*CLHEP::barn , 39.95*CLHEP::barn , 350 91.15*CLHEP::barn , 104.4*CLHEP::barn , 351 352 static const G4double Tdat[22] = { 353 100*CLHEP::eV, 200*CLHEP::eV, 400*CLHEP: 354 1*CLHEP::keV, 2*CLHEP::keV, 4*CLHEP::k 355 10*CLHEP::keV, 20*CLHEP::keV, 40*CLHEP:: 356 100*CLHEP::keV, 200*CLHEP::keV, 400*CLHEP: 357 1*CLHEP::MeV, 2*CLHEP::MeV, 4*CLHEP::M 358 10*CLHEP::MeV, 20*CLHEP::MeV}; 359 360 if(eKineticEnergy <= Tlim) 361 { 362 // get bin number in T (beta2) 363 G4int iT = 21; 364 // Loop checking, 03-Aug-2015, Vladimir Iv 365 while ((iT>=0)&&(Tdat[iT]>=eKineticEnergy) 366 367 iT = std::min(std::max(iT, 0), 20); 368 369 // calculate betasquare values 370 G4double T = Tdat[iT]; 371 G4double E = T + CLHEP::electron_mass_c2; 372 G4double b2small = T*(E+CLHEP::electron_ma 373 374 T = Tdat[iT+1]; 375 E = T + CLHEP::electron_mass_c2; 376 G4double b2big = T*(E+CLHEP::electron_mass 377 G4double ratb2 = (beta2-b2small)/(b2big-b2 378 379 if (charge < 0.) 380 { 381 c1 = celectron[iZ][iT]; 382 c2 = celectron[iZ+1][iT]; 383 cc1 = c1+ratZ*(c2-c1); 384 385 c1 = celectron[iZ][iT+1]; 386 c2 = celectron[iZ+1][iT+1]; 387 } 388 else 389 { 390 c1 = cpositron[iZ][iT]; 391 c2 = cpositron[iZ+1][iT]; 392 cc1 = c1+ratZ*(c2-c1); 393 394 c1 = cpositron[iZ][iT+1]; 395 c2 = cpositron[iZ+1][iT+1]; 396 } 397 G4double cc2 = c1+ratZ*(c2-c1); 398 sigma *= sigmafactor/(cc1+ratb2*(cc2-cc1)) 399 } 400 else 401 { 402 c1 = bg2lim*sig0[iZ]*(1.+hecorr[iZ]*(beta2 403 c2 = bg2lim*sig0[iZ+1]*(1.+hecorr[iZ+1]*(b 404 if((atomicNumber >= ZZ1) && (atomicNumber 405 sigma = c1+ratZ*(c2-c1) ; 406 else if(atomicNumber < ZZ1) 407 sigma = atomicNumber*atomicNumber*c1/(ZZ 408 else if(atomicNumber > ZZ2) 409 sigma = atomicNumber*atomicNumber*c2/(ZZ 410 } 411 // low energy correction based on theory 412 sigma *= (1.+0.30/(1.+std::sqrt(1000.*eKinet 413 414 return sigma; 415 } 416 417 //....oooOO0OOooo........oooOO0OOooo........oo 418 419 void G4UrbanMscModel::StartTracking(G4Track* t 420 { 421 SetParticle(track->GetDynamicParticle()->Get 422 firstStep = true; 423 insideskin = false; 424 fr = facrange; 425 tlimit = tgeom = rangeinit = geombig; 426 smallstep = 1.e10; 427 stepmin = tlimitminfix; 428 tlimitmin = 10.*tlimitminfix; 429 rndmEngineMod = G4Random::getTheEngine(); 430 } 431 432 //....oooOO0OOooo........oooOO0OOooo........oo 433 434 G4double G4UrbanMscModel::ComputeTruePathLengt 435 const G4Track& tr 436 G4double& current 437 { 438 tPathLength = currentMinimalStep; 439 const G4DynamicParticle* dp = track.GetDynam 440 441 G4StepPoint* sp = track.GetStep()->GetPreSte 442 G4StepStatus stepStatus = sp->GetStepStatus( 443 couple = track.GetMaterialCutsCouple(); 444 SetCurrentCouple(couple); 445 idx = couple->GetIndex(); 446 currentKinEnergy = dp->GetKineticEnergy(); 447 currentLogKinEnergy = dp->GetLogKineticEnerg 448 currentRange = GetRange(particle,currentKinE 449 lambda0 = GetTransportMeanFreePath(particle, 450 451 tPathLength = std::min(tPathLength,currentRa 452 /* 453 G4cout << "G4Urban::StepLimit tPathLength= " 454 << " range= " <<currentRange<< " lambda= "<< 455 <<G4endl; 456 */ 457 // extreme small step 458 if(tPathLength < tlimitminfix) { 459 latDisplasment = false; 460 return ConvertTrueToGeom(tPathLength, curr 461 } 462 463 presafety = (stepStatus == fGeomBoundary) ? 464 : ComputeSafety(sp->GetPosition( 465 466 // stop here if small step or range is less 467 if((tPathLength == currentRange && tPathLeng 468 tPathLength < tlimitminfix) { 469 latDisplasment = false; 470 return ConvertTrueToGeom(tPathLength, curr 471 } 472 473 // upper limit for the straight line distanc 474 // for electrons and positrons 475 G4double distance = (mass < masslimite) 476 ? currentRange*msc[idx]->doverra 477 // for muons, hadrons 478 : currentRange*msc[idx]->doverrb; 479 480 /* 481 G4cout << "G4Urban::StepLimit tPathLength= " 482 <<tPathLength<<" safety= " << pres 483 << " range= " <<currentRange<< " lam 484 << " Alg: " << steppingAlgorithm < 485 */ 486 // far from geometry boundary 487 if(distance < presafety) 488 { 489 latDisplasment = false; 490 return ConvertTrueToGeom(tPathLength, cu 491 } 492 493 latDisplasment = latDisplasmentbackup; 494 // ----------------------------------------- 495 // distance to boundary 496 if (steppingAlgorithm == fUseDistanceToBound 497 { 498 //compute geomlimit and presafety 499 geomlimit = ComputeGeomLimit(track, pres 500 /* 501 G4cout << "G4Urban::Distance to bounda 502 <<geomlimit<<" safety= " << presaf 503 */ 504 505 smallstep += 1.; 506 insideskin = false; 507 tgeom = geombig; 508 509 // initialisation at first step and at t 510 if(firstStep || (stepStatus == fGeomBoun 511 { 512 rangeinit = currentRange; 513 if(!firstStep) { smallstep = 1.; } 514 515 //stepmin ~ lambda_elastic 516 stepmin = ComputeStepmin(); 517 skindepth = skin*stepmin; 518 tlimitmin = ComputeTlimitmin(); 519 /* 520 G4cout << "rangeinit= " << rangeinit 521 << " tlimitmin= " << tlimitmi 522 << geomlimit <<G4endl; 523 */ 524 } 525 // constraint from the geometry 526 if((geomlimit < geombig) && (geomlimit > 527 { 528 // geomlimit is a geometrical step l 529 // transform it to true path length 530 if(lambda0 > geomlimit) { 531 geomlimit = -lambda0*G4Log(1.-geom 532 } 533 tgeom = (stepStatus == fGeomBoundary 534 : facrange*rangeinit + stepmin; 535 } 536 537 //step limit 538 tlimit = (currentRange > presafety) ? 539 std::max(facrange*rangeinit, facsafety 540 541 //lower limit for tlimit 542 tlimit = std::min(std::max(tlimit,tlimit 543 /* 544 G4cout << "tgeom= " << tgeom << " geomli 545 << " tlimit= " << tlimit << " pres 546 */ 547 // shortcut 548 if((tPathLength < tlimit) && (tPathLengt 549 (smallstep > skin) && (tPathLength < 550 { 551 return ConvertTrueToGeom(tPathLength, 552 } 553 554 // step reduction near to boundary 555 if(smallstep <= skin) 556 { 557 tlimit = stepmin; 558 insideskin = true; 559 } 560 else if(geomlimit < geombig) 561 { 562 if(geomlimit > skindepth) 563 { 564 tlimit = std::min(tlimit, geomli 565 } 566 else 567 { 568 insideskin = true; 569 tlimit = std::min(tlimit, stepmi 570 } 571 } 572 573 tlimit = std::max(tlimit, stepmin); 574 575 // randomise if not 'small' step and ste 576 tPathLength = (tlimit < tPathLength && s 577 ? std::min(tPathLength, Randomizetlimi 578 : std::min(tPathLength, tlimit); 579 } 580 // ----------------------------------------- 581 // for simulation with or without magnetic f 582 // there no small step/single scattering at 583 else if(steppingAlgorithm == fUseSafety) 584 { 585 if(firstStep || (stepStatus == fGeomBoun 586 rangeinit = currentRange; 587 fr = facrange; 588 // stepping for e+/e- only (not for mu 589 if(mass < masslimite) 590 { 591 rangeinit = std::max(rangeinit, la 592 if(lambda0 > lambdalimit) { 593 fr *= (0.75+0.25*lambda0/lambdal 594 } 595 } 596 //lower limit for tlimit 597 stepmin = ComputeStepmin(); 598 tlimitmin = ComputeTlimitmin(); 599 } 600 601 //step limit 602 tlimit = (currentRange > presafety) ? 603 std::max(fr*rangeinit, facsafety*presa 604 605 //lower limit for tlimit 606 tlimit = std::max(tlimit, tlimitmin); 607 608 // randomise if step determined by msc 609 tPathLength = (tlimit < tPathLength) ? 610 std::min(tPathLength, Randomizetlimit( 611 } 612 // ----------------------------------------- 613 // for simulation with or without magnetic f 614 // there is small step/single scattering at 615 else if(steppingAlgorithm == fUseSafetyPlus) 616 { 617 if(firstStep || (stepStatus == fGeomBoun 618 rangeinit = currentRange; 619 fr = facrange; 620 if(mass < masslimite) 621 { 622 if(lambda0 > lambdalimit) { 623 fr *= (0.84+0.16*lambda0/lambdal 624 } 625 } 626 //lower limit for tlimit 627 stepmin = ComputeStepmin(); 628 tlimitmin = ComputeTlimitmin(); 629 } 630 //step limit 631 tlimit = (currentRange > presafety) ? 632 std::max(fr*rangeinit, facsafety*presafety) 633 634 //lower limit for tlimit 635 tlimit = std::max(tlimit, tlimitmin); 636 637 // condition for tPathLength from drr an 638 if(currentRange > finalr) { 639 G4double tmax = drr*currentRange+ 640 finalr*(1.-drr)*(2.-fi 641 tPathLength = std::min(tPathLength,tma 642 } 643 644 // randomise if step determined by msc 645 tPathLength = (tlimit < tPathLength) ? 646 std::min(tPathLength, Randomizetlimit( 647 } 648 649 // ----------------------------------------- 650 // simple step limitation 651 else 652 { 653 if (stepStatus == fGeomBoundary) 654 { 655 tlimit = (currentRange > lambda0) 656 ? facrange*currentRange : facrange*lambd 657 tlimit = std::max(tlimit, tlimitmin) 658 } 659 // randomise if step determined by msc 660 tPathLength = (tlimit < tPathLength) ? 661 std::min(tPathLength, Randomizetlimit( 662 } 663 664 // ----------------------------------------- 665 firstStep = false; 666 return ConvertTrueToGeom(tPathLength, curren 667 } 668 669 //....oooOO0OOooo........oooOO0OOooo........oo 670 671 G4double G4UrbanMscModel::ComputeGeomPathLengt 672 { 673 lambdaeff = lambda0; 674 par1 = -1. ; 675 par2 = par3 = 0. ; 676 677 // this correction needed to run MSC with eI 678 // and makes no harm for a normal run 679 tPathLength = std::min(tPathLength,currentRa 680 681 // do the true -> geom transformation 682 zPathLength = tPathLength; 683 684 // z = t for very small tPathLength 685 if(tPathLength < tlimitminfix2) return zPath 686 687 /* 688 G4cout << "ComputeGeomPathLength: tpl= " << 689 << " R= " << currentRange << " L0= " 690 << " E= " << currentKinEnergy << " " 691 << particle->GetParticleName() << G4e 692 */ 693 G4double tau = tPathLength/lambda0 ; 694 695 if ((tau <= tausmall) || insideskin) { 696 zPathLength = std::min(tPathLength, lambda 697 698 } else if (tPathLength < currentRange*dtrl) 699 zPathLength = (tau < taulim) ? tPathLength 700 : lambda0*(1.-G4Exp(-tau)); 701 702 } else if(currentKinEnergy < mass || tPathLe 703 par1 = 1./currentRange; 704 par2 = currentRange/lambda0; 705 par3 = 1.+par2; 706 if(tPathLength < currentRange) { 707 zPathLength = 708 (1.-G4Exp(par3*G4Log(1.-tPathLength/cu 709 } else { 710 zPathLength = 1./(par1*par3); 711 } 712 713 } else { 714 G4double rfin = std::max(currentRange-tPat 715 G4double T1 = GetEnergy(particle,rfin,coup 716 G4double lambda1 = GetTransportMeanFreePat 717 718 par1 = (lambda0-lambda1)/(lambda0*tPathLen 719 //G4cout << "par1= " << par1 << " L1= " << 720 par2 = 1./(par1*lambda0); 721 par3 = 1.+par2; 722 zPathLength = (1.-G4Exp(par3*G4Log(lambda1 723 } 724 725 zPathLength = std::min(zPathLength, lambda0) 726 //G4cout<< "zPathLength= "<< zPathLength<< " 727 return zPathLength; 728 } 729 730 //....oooOO0OOooo........oooOO0OOooo........oo 731 732 G4double G4UrbanMscModel::ComputeTrueStepLengt 733 { 734 // step defined other than transportation 735 if(geomStepLength == zPathLength) { 736 //G4cout << "Urban::ComputeTrueLength: tPa 737 // << " step= " << geomStepLengt 738 return tPathLength; 739 } 740 741 zPathLength = geomStepLength; 742 743 // t = z for very small step 744 if(geomStepLength < tlimitminfix2) { 745 tPathLength = geomStepLength; 746 747 // recalculation 748 } else { 749 750 G4double tlength = geomStepLength; 751 if((geomStepLength > lambda0*tausmall) && 752 753 if(par1 < 0.) { 754 tlength = -lambda0*G4Log(1.-geomStepLe 755 } else { 756 const G4double par4 = par1*par3; 757 if(par4*geomStepLength < 1.) { 758 tlength = (1.-G4Exp(G4Log(1.-par4*ge 759 } else { 760 tlength = currentRange; 761 } 762 } 763 764 if(tlength < geomStepLength) { tlength 765 else if(tlength > tPathLength) { tlength 766 } 767 tPathLength = tlength; 768 } 769 //G4cout << "Urban::ComputeTrueLength: tPath 770 // << " step= " << geomStepLength << 771 772 return tPathLength; 773 } 774 775 //....oooOO0OOooo........oooOO0OOooo........oo 776 777 G4ThreeVector& 778 G4UrbanMscModel::SampleScattering(const G4Thre 779 G4double /*s 780 { 781 fDisplacement.set(0.0,0.0,0.0); 782 if(tPathLength >= currentRange) { return fDi 783 784 G4double kinEnergy = currentKinEnergy; 785 if (tPathLength > currentRange*dtrl) { 786 kinEnergy = GetEnergy(particle,currentRang 787 } else if(tPathLength > currentRange*0.01) { 788 kinEnergy -= tPathLength*GetDEDX(particle, 789 currentLo 790 } 791 792 if((tPathLength <= tlimitminfix) || (tPathLe 793 (kinEnergy <= CLHEP::eV)) { return fDispl 794 795 G4double cth = SampleCosineTheta(tPathLength 796 797 // protection against 'bad' cth values 798 if(std::abs(cth) >= 1.0) { return fDisplacem 799 800 G4double sth = std::sqrt((1.0 - cth)*(1.0 + 801 G4double phi = CLHEP::twopi*rndmEngineMod->f 802 G4ThreeVector newDirection(sth*std::cos(phi) 803 newDirection.rotateUz(oldDirection); 804 805 fParticleChange->ProposeMomentumDirection(ne 806 /* 807 G4cout << "G4UrbanMscModel::SampleSecondarie 808 << " sinTheta= " << sth << " safety(m 809 << " trueStep(mm)= " << tPathLength 810 << " geomStep(mm)= " << zPathLength 811 << G4endl; 812 */ 813 814 if (latDisplasment && currentTau >= tausmall 815 if(dispAlg96) { SampleDisplacement(sth, ph 816 else { SampleDisplacementNew(cth, 817 fDisplacement.rotateUz(oldDirection); 818 } 819 return fDisplacement; 820 } 821 822 //....oooOO0OOooo........oooOO0OOooo........oo 823 824 G4double G4UrbanMscModel::SampleCosineTheta(G4 825 G4 826 { 827 G4double cth = 1.0; 828 G4double tau = trueStepLength/lambda0; 829 830 // mean tau value 831 if(currentKinEnergy != kinEnergy) { 832 G4double lambda1 = GetTransportMeanFreePat 833 if(std::abs(lambda1 - lambda0) > lambda0*0 834 tau = trueStepLength*G4Log(lambda0/lambd 835 } 836 } 837 838 currentTau = tau; 839 lambdaeff = trueStepLength/currentTau; 840 currentRadLength = couple->GetMaterial()->Ge 841 842 if (tau >= taubig) { cth = -1.+2.*rndmEngine 843 else if (tau >= tausmall) { 844 static const G4double numlim = 0.01; 845 static const G4double onethird = 1./3.; 846 if(tau < numlim) { 847 xmeanth = 1.0 - tau*(1.0 - 0.5*tau); 848 x2meanth= 1.0 - tau*(5.0 - 6.25*tau)*one 849 } else { 850 xmeanth = G4Exp(-tau); 851 x2meanth = (1.+2.*G4Exp(-2.5*tau))*oneth 852 } 853 854 // too large step of low-energy particle 855 G4double relloss = 1. - kinEnergy/currentK 856 static const G4double rellossmax= 0.50; 857 if(relloss > rellossmax) { 858 return SimpleScattering(); 859 } 860 // is step extreme small ? 861 G4bool extremesmallstep = false; 862 G4double tsmall = std::min(tlimitmin,lambd 863 864 G4double theta0; 865 if(trueStepLength > tsmall) { 866 theta0 = ComputeTheta0(trueStepLength,ki 867 } else { 868 theta0 = std::sqrt(trueStepLength/tsmall 869 *ComputeTheta0(tsmall,kinEnergy); 870 extremesmallstep = true; 871 } 872 873 static const G4double onesixth = 1./6.; 874 static const G4double one12th = 1./12.; 875 static const G4double theta0max = CLHEP::p 876 877 // protection for very small angles 878 G4double theta2 = theta0*theta0; 879 880 if(theta2 < tausmall) { return cth; } 881 if(theta0 > theta0max) { return SimpleScat 882 883 G4double x = theta2*(1.0 - theta2*one12th) 884 if(theta2 > numlim) { 885 G4double sth = 2*std::sin(0.5*theta0); 886 x = sth*sth; 887 } 888 889 // parameter for tail 890 G4double ltau = G4Log(tau); 891 G4double u = !extremesmallstep ? G4Exp(lta 892 : G4Exp(G4Log(tsmall/lambda0)*onesixth); 893 894 G4double xx = G4Log(lambdaeff/currentRadL 895 G4double xsi = msc[idx]->coeffc1 + 896 u*(msc[idx]->coeffc2+msc[idx]->coeffc3*u 897 898 // tail should not be too big 899 xsi = std::max(xsi, 1.9); 900 /* 901 if(KineticEnergy > 20*MeV && xsi < 1.6) 902 G4cout << "G4UrbanMscModel::SampleCosi 903 << KineticEnergy/GeV 904 << " !!** c= " << xsi 905 << " **!! length(mm)= " << true 906 << " " << couple->GetMaterial() 907 << " tau= " << tau << G4endl; 908 } 909 */ 910 911 G4double c = xsi; 912 913 if(std::abs(c-3.) < 0.001) { c = 3.00 914 else if(std::abs(c-2.) < 0.001) { c = 2.00 915 916 G4double c1 = c-1.; 917 G4double ea = G4Exp(-xsi); 918 G4double eaa = 1.-ea ; 919 G4double xmean1 = 1.-(1.-(1.+xsi)*ea)*x/ea 920 G4double x0 = 1. - xsi*x; 921 922 // G4cout << " xmean1= " << xmean1 << " x 923 924 if(xmean1 <= 0.999*xmeanth) { return Simpl 925 926 //from continuity of derivatives 927 G4double b = 1.+(c-xsi)*x; 928 929 G4double b1 = b+1.; 930 G4double bx = c*x; 931 932 G4double eb1 = G4Exp(G4Log(b1)*c1); 933 G4double ebx = G4Exp(G4Log(bx)*c1); 934 G4double d = ebx/eb1; 935 936 G4double xmean2 = (x0 + d - (bx - b1*d)/(c 937 938 G4double f1x0 = ea/eaa; 939 G4double f2x0 = c1/(c*(1. - d)); 940 G4double prob = f2x0/(f1x0+f2x0); 941 942 G4double qprob = xmeanth/(prob*xmean1+(1.- 943 944 // sampling of costheta 945 //G4cout << "c= " << c << " qprob= " << qp 946 // << " c1= " << c1 << " b1= " << b1 << " 947 // << G4endl; 948 rndmEngineMod->flatArray(2, rndmarray); 949 if(rndmarray[0] < qprob) 950 { 951 G4double var = 0; 952 if(rndmarray[1] < prob) { 953 cth = 1.+G4Log(ea+rndmEngineMod->flat( 954 } else { 955 var = (1.0 - d)*rndmEngineMod->flat(); 956 if(var < numlim*d) { 957 var /= (d*c1); 958 cth = -1.0 + var*(1.0 - 0.5*var*c)*( 959 } else { 960 cth = 1. + x*(c - xsi - c*G4Exp(-G4L 961 } 962 } 963 } else { 964 cth = -1.+2.*rndmarray[1]; 965 } 966 } 967 return cth; 968 } 969 970 //....oooOO0OOooo........oooOO0OOooo........oo 971 972 G4double G4UrbanMscModel::ComputeTheta0(G4doub 973 G4doub 974 { 975 // for all particles take the width of the c 976 // from a parametrization similar to the H 977 // ( Highland formula: Particle Physics Book 978 G4double invbetacp = (kinEnergy+mass)/(kinEn 979 if(currentKinEnergy != kinEnergy) { 980 invbetacp = std::sqrt(invbetacp*(currentKi 981 (currentKinEnergy*(currentKinEnergy+2. 982 } 983 G4double y = trueStepLength/currentRadLength 984 985 if(fPosiCorrection && particle == positron) 986 { 987 static const G4double xl= 0.6; 988 static const G4double xh= 0.9; 989 static const G4double e = 113.0; 990 G4double corr; 991 992 G4double tau = std::sqrt(currentKinEnergy* 993 G4double x = std::sqrt(tau*(tau+2.)/((tau+ 994 G4double a = msc[idx]->posa; 995 G4double b = msc[idx]->posb; 996 G4double c = msc[idx]->posc; 997 G4double d = msc[idx]->posd; 998 if(x < xl) { 999 corr = a*(1.-G4Exp(-b*x)); 1000 } else if(x > xh) { 1001 corr = c+d*G4Exp(e*(x-1.)); 1002 } else { 1003 G4double yl = a*(1.-G4Exp(-b*xl)); 1004 G4double yh = c+d*G4Exp(e*(xh-1.)); 1005 G4double y0 = (yh-yl)/(xh-xl); 1006 G4double y1 = yl-y0*xl; 1007 corr = y0*x+y1; 1008 } 1009 //======================================= 1010 y *= corr*msc[idx]->pose; 1011 } 1012 1013 static const G4double c_highland = 13.6*CLH 1014 G4double theta0 = c_highland*std::abs(charg 1015 1016 // correction factor from e- scattering dat 1017 theta0 *= (msc[idx]->coeffth1+msc[idx]->coe 1018 return theta0; 1019 } 1020 1021 //....oooOO0OOooo........oooOO0OOooo........o 1022 1023 void G4UrbanMscModel::SampleDisplacement(G4do 1024 { 1025 // simple and fast sampling 1026 // based on single scattering results 1027 // u = r/rmax : mean value 1028 1029 G4double rmax = std::sqrt((tPathLength-zPat 1030 if(rmax > 0.) 1031 { 1032 G4double r = 0.73*rmax; 1033 1034 // simple distribution for v=Phi-phi=psi 1035 // beta determined from the requirement t 1036 // the same mean value than that obtained 1037 1038 static const G4double cbeta = 2.160; 1039 static const G4double cbeta1 = 1. - G4Exp 1040 rndmEngineMod->flatArray(2, rndmarray); 1041 G4double psi = -G4Log(1. - rndmarray[0]*c 1042 G4double Phi = (rndmarray[1] < 0.5) ? phi 1043 fDisplacement.set(r*std::cos(Phi),r*std:: 1044 } 1045 } 1046 1047 //....oooOO0OOooo........oooOO0OOooo........o 1048 1049 void G4UrbanMscModel::SampleDisplacementNew(G 1050 { 1051 // simple and fast sampling 1052 // based on single scattering results 1053 // u = (r/rmax)**2 : distribution from ss s 1054 const G4double eps = 1.e-3; 1055 const G4double rmax = 1056 std::sqrt((tPathLength-zPathLength)*(tPat 1057 1058 if(rmax > 0.) 1059 { 1060 const G4double x0 = 0.73 ; 1061 const G4double alpha = G4Log(7.33)/x0 ; 1062 const G4double a1 = 1.-x0 ; 1063 const G4double a2 = 1.-G4Exp(-alpha*x0) ; 1064 const G4double a3 = G4Exp(alpha*x0)-1. ; 1065 const G4double w1 = 2.*a2/(alpha*a1+2.*a2 1066 1067 G4double r, sqx; 1068 if (rmax/currentRange < eps) 1069 { 1070 r = 0.73*rmax ; 1071 sqx = 1.; 1072 } 1073 else 1074 { 1075 rndmEngineMod->flatArray(2,rndmarray); 1076 const G4double x = (rndmarray[0] < w1) 1077 1. - a1*std::sqrt(1.-rndmarray[1]); 1078 1079 sqx = std::sqrt(x); 1080 r = sqx*rmax; 1081 } 1082 // Gaussian distribution for Phi-phi=psi 1083 const G4double sigma = 0.1+0.9*sqx; 1084 const G4double psi = G4RandGauss::shoot(0 1085 const G4double Phi = phi+psi; 1086 fDisplacement.set(r*std::cos(Phi), r*std: 1087 } 1088 } 1089 1090 //....oooOO0OOooo........oooOO0OOooo........o 1091 1092 void G4UrbanMscModel::InitialiseModelCache() 1093 { 1094 // it is assumed, that for the second run o 1095 // of a new G4MaterialCutsCouple is possibl 1096 auto theCoupleTable = G4ProductionCutsTable 1097 std::size_t numOfCouples = theCoupleTable-> 1098 if(numOfCouples != msc.size()) { msc.resize 1099 1100 for(G4int j=0; j<(G4int)numOfCouples; ++j) 1101 auto aCouple = theCoupleTable->GetMateria 1102 1103 // new couple 1104 msc[j] = new mscData(); 1105 G4double Zeff = aCouple->GetMaterial()->G 1106 G4double sqrz = std::sqrt(Zeff); 1107 msc[j]->sqrtZ = sqrz; 1108 // parameterisation of step limitation 1109 msc[j]->factmin = dispAlg96 ? 0.001 : 0.0 1110 G4double lnZ = G4Log(Zeff); 1111 // correction in theta0 formula 1112 G4double w = G4Exp(lnZ/6.); 1113 G4double facz = 0.990395+w*(-0.168386+w*0 1114 msc[j]->coeffth1 = facz*(1. - 8.7780e-2/Z 1115 msc[j]->coeffth2 = facz*(4.0780e-2 + 1.73 1116 1117 // tail parameters 1118 G4double Z13 = w*w; 1119 msc[j]->coeffc1 = 2.3785 - Z13*(4.1981 1120 msc[j]->coeffc2 = 4.7526e-1 + Z13*(1.7694 1121 msc[j]->coeffc3 = 2.3683e-1 - Z13*(1.8111 1122 msc[j]->coeffc4 = 1.7888e-2 + Z13*(1.9659 1123 1124 msc[j]->Z23 = Z13*Z13; 1125 1126 msc[j]->stepmina = 27.725/(1.+0.203*Zeff) 1127 msc[j]->stepminb = 6.152/(1.+0.111*Zeff) 1128 1129 // 21.07.2020 1130 msc[j]->doverra = 9.6280e-1 - 8.4848e-2*m 1131 1132 // 06.10.2020 1133 // msc[j]->doverra = 7.7024e-1 - 6.7878e- 1134 msc[j]->doverrb = 1.15 - 9.76e-4*Zeff; 1135 1136 // corrections for e+ 1137 msc[j]->posa = 0.994-4.08e-3*Zeff; 1138 msc[j]->posb = 7.16+(52.6+365./Zeff)/Zeff 1139 msc[j]->posc = 1.000-4.47e-3*Zeff; 1140 msc[j]->posd = 1.21e-3*Zeff; 1141 msc[j]->pose = 1.+Zeff*(1.84035e-4*Zeff-1 1142 } 1143 } 1144 1145 //....oooOO0OOooo........oooOO0OOooo........o 1146