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 // GEANT 4 class implementation file 28 // 29 // History: first implementation, 30 // 21-5-98 V.Grichine 31 // 28-05-01, V.Ivanchenko minor changes t 32 // 04.03.05, V.Grichine: get local field 33 // 19-05-06, V.Ivanchenko rename from G4S 34 // 35 ////////////////////////////////////////////// 36 37 #include "G4SynchrotronRadiationInMat.hh" 38 39 #include "G4EmProcessSubType.hh" 40 #include "G4Field.hh" 41 #include "G4FieldManager.hh" 42 #include "G4Integrator.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4PropagatorInField.hh" 45 #include "G4SystemOfUnits.hh" 46 #include "G4PhysicsModelCatalog.hh" 47 48 const G4double G4SynchrotronRadiationInMat::fI 49 1.000000e+00, 9.428859e-01, 9.094095e-01, 8. 50 8.337008e-01, 8.124961e-01, 7.925217e-01, 7. 51 7.381233e-01, 7.214521e-01, 7.053634e-01, 6. 52 6.600922e-01, 6.458793e-01, 6.320533e-01, 6. 53 5.926459e-01, 5.801347e-01, 5.679103e-01, 5. 54 5.328395e-01, 5.216482e-01, 5.106904e-01, 4. 55 4.791351e-01, 4.690316e-01, 4.591249e-01, 4. 56 4.305320e-01, 4.213608e-01, 4.123623e-01, 4. 57 3.863639e-01, 3.780179e-01, 3.698262e-01, 3. 58 3.461460e-01, 3.385411e-01, 3.310757e-01, 3. 59 3.094921e-01, 3.025605e-01, 2.957566e-01, 2. 60 2.760907e-01, 2.697773e-01, 2.635817e-01, 2. 61 2.456834e-01, 2.399409e-01, 2.343074e-01, 2. 62 2.180442e-01, 2.128303e-01, 2.077174e-01, 2. 63 1.929696e-01, 1.882457e-01, 1.836155e-01, 1. 64 1.702730e-01, 1.660036e-01, 1.618212e-01, 1. 65 1.497822e-01, 1.459344e-01, 1.421671e-01, 1. 66 1.313360e-01, 1.278785e-01, 1.244956e-01, 1. 67 1.147818e-01, 1.116850e-01, 1.086570e-01, 1. 68 9.997405e-02, 9.720975e-02, 9.450865e-02, 9. 69 8.677391e-02, 8.431501e-02, 8.191406e-02, 7. 70 7.504872e-02, 7.286944e-02, 7.074311e-02, 6. 71 6.467208e-02, 6.274790e-02, 6.087191e-02, 5. 72 5.552387e-02, 5.383150e-02, 5.218282e-02, 5. 73 4.749020e-02, 4.600763e-02, 4.456450e-02, 4. 74 4.046353e-02, 3.917002e-02, 3.791195e-02, 3. 75 3.434274e-02, 3.321884e-02, 3.212665e-02, 3. 76 2.903319e-02, 2.806076e-02, 2.711656e-02, 2. 77 2.444677e-02, 2.360897e-02, 2.279620e-02, 2. 78 2.050194e-02, 1.978324e-02, 1.908662e-02, 1. 79 1.712363e-02, 1.650979e-02, 1.591533e-02, 1. 80 1.424314e-02, 1.372117e-02, 1.321613e-02, 1. 81 1.179798e-02, 1.135611e-02, 1.092896e-02, 1. 82 9.731635e-03, 9.359254e-03, 8.999595e-03, 8. 83 7.993280e-03, 7.680879e-03, 7.379426e-03, 7. 84 6.537491e-03, 6.276605e-03, 6.025092e-03, 5. 85 5.323912e-03, 5.107045e-03, 4.898164e-03, 4. 86 4.316896e-03, 4.137454e-03, 3.964780e-03, 3. 87 3.485150e-03, 3.337364e-03, 3.195284e-03, 3. 88 2.801361e-03, 2.680213e-03, 2.563852e-03, 2. 89 }; 90 91 ////////////////////////////////////////////// 92 // Constructor 93 G4SynchrotronRadiationInMat::G4SynchrotronRadi 94 const G4String& processName, G4ProcessType t 95 : G4VDiscreteProcess(processName, type) 96 , theGamma(G4Gamma::Gamma()) 97 , theElectron(G4Electron::Electron()) 98 , thePositron(G4Positron::Positron()) 99 , LowestKineticEnergy(10. * keV) 100 , fAlpha(0.0) 101 , fRootNumber(80) 102 , fVerboseLevel(verboseLevel) 103 { 104 G4TransportationManager* transportMgr = 105 G4TransportationManager::GetTransportation 106 107 fFieldPropagator = transportMgr->GetPropagat 108 secID = G4PhysicsModelCatalog::GetModelID("m 109 SetProcessSubType(fSynchrotronRadiation); 110 CutInRange = GammaCutInKineticEnergyNow = El 111 PositronCutInKineticEnergyNow = ParticleCu 112 fPsiGamma = fEta = fOrderAngleK = 0.0; 113 } 114 115 ////////////////////////////////////////////// 116 // Destructor 117 G4SynchrotronRadiationInMat::~G4SynchrotronRad 118 119 G4bool G4SynchrotronRadiationInMat::IsApplicab 120 const G4ParticleDefinition& particle) 121 { 122 return ((&particle == (const G4ParticleDefin 123 (&particle == (const G4ParticleDefin 124 } 125 126 G4double G4SynchrotronRadiationInMat::GetLambd 127 128 G4double G4SynchrotronRadiationInMat::GetEnerg 129 130 // Production of synchrotron X-ray photon 131 // Geant4 internal units. 132 G4double G4SynchrotronRadiationInMat::GetMeanF 133 const G4Track& trackData, G4double, G4ForceC 134 { 135 // gives the MeanFreePath in GEANT4 internal 136 G4double MeanFreePath; 137 138 const G4DynamicParticle* aDynamicParticle = 139 140 *condition = NotForced; 141 142 G4double gamma = 143 aDynamicParticle->GetTotalEnergy() / aDyna 144 145 G4double particleCharge = aDynamicParticle-> 146 147 G4double KineticEnergy = aDynamicParticle->G 148 149 if(KineticEnergy < LowestKineticEnergy || ga 150 MeanFreePath = DBL_MAX; 151 else 152 { 153 G4ThreeVector FieldValue; 154 const G4Field* pField = nullptr; 155 156 G4FieldManager* fieldMgr = nullptr; 157 G4bool fieldExertsForce = false; 158 159 if((particleCharge != 0.0)) 160 { 161 fieldMgr = 162 fFieldPropagator->FindAndSetFieldManag 163 164 if(fieldMgr != nullptr) 165 { 166 // If the field manager has no field, 167 fieldExertsForce = (fieldMgr->GetDetec 168 } 169 } 170 if(fieldExertsForce) 171 { 172 pField = fieldMgr->G 173 G4ThreeVector globPosition = trackData.G 174 175 G4double globPosVec[4], FieldValueVec[6] 176 177 globPosVec[0] = globPosition.x(); 178 globPosVec[1] = globPosition.y(); 179 globPosVec[2] = globPosition.z(); 180 globPosVec[3] = trackData.GetGlobalTime( 181 182 pField->GetFieldValue(globPosVec, FieldV 183 184 FieldValue = 185 G4ThreeVector(FieldValueVec[0], FieldV 186 187 G4ThreeVector unitMomentum = aDynamicPar 188 G4ThreeVector unitMcrossB = FieldValue. 189 G4double perpB = unitMcrossB 190 G4double beta = aDynamicPar 191 (aDynamicParticle->GetTo 192 193 if(perpB > 0.0) 194 MeanFreePath = fLambdaConst * beta / p 195 else 196 MeanFreePath = DBL_MAX; 197 } 198 else 199 MeanFreePath = DBL_MAX; 200 } 201 if(fVerboseLevel > 0) 202 { 203 G4cout << "G4SynchrotronRadiationInMat::Me 204 << " m" << G4endl; 205 } 206 return MeanFreePath; 207 } 208 209 ////////////////////////////////////////////// 210 G4VParticleChange* G4SynchrotronRadiationInMat 211 const G4Track& trackData, const G4Step& step 212 213 { 214 aParticleChange.Initialize(trackData); 215 216 const G4DynamicParticle* aDynamicParticle = 217 218 G4double gamma = 219 aDynamicParticle->GetTotalEnergy() / (aDyn 220 221 if(gamma <= 1.0e3) 222 { 223 return G4VDiscreteProcess::PostStepDoIt(tr 224 } 225 G4double particleCharge = aDynamicParticle-> 226 227 G4ThreeVector FieldValue; 228 const G4Field* pField = nullptr; 229 230 G4FieldManager* fieldMgr = nullptr; 231 G4bool fieldExertsForce = false; 232 233 if((particleCharge != 0.0)) 234 { 235 fieldMgr = fFieldPropagator->FindAndSetFie 236 if(fieldMgr != nullptr) 237 { 238 // If the field manager has no field, th 239 fieldExertsForce = (fieldMgr->GetDetecto 240 } 241 } 242 if(fieldExertsForce) 243 { 244 pField = fieldMgr->Get 245 G4ThreeVector globPosition = trackData.Get 246 G4double globPosVec[4], FieldValueVec[6]; 247 globPosVec[0] = globPosition.x(); 248 globPosVec[1] = globPosition.y(); 249 globPosVec[2] = globPosition.z(); 250 globPosVec[3] = trackData.GetGlobalTime(); 251 252 pField->GetFieldValue(globPosVec, FieldVal 253 FieldValue = 254 G4ThreeVector(FieldValueVec[0], FieldVal 255 256 G4ThreeVector unitMomentum = aDynamicParti 257 G4ThreeVector unitMcrossB = FieldValue.cr 258 G4double perpB = unitMcrossB.m 259 if(perpB > 0.0) 260 { 261 // M-C of synchrotron photon energy 262 G4double energyOfSR = GetRandomEnergySR( 263 264 if(fVerboseLevel > 0) 265 { 266 G4cout << "SR photon energy = " << ene 267 } 268 269 // check against insufficient energy 270 if(energyOfSR <= 0.0) 271 { 272 return G4VDiscreteProcess::PostStepDoI 273 } 274 G4double kineticEnergy = aDynamicParticl 275 G4ParticleMomentum particleDirection = 276 aDynamicParticle->GetMomentumDirection 277 278 // M-C of its direction, simplified dipo 279 G4double cosTheta, sinTheta, fcos, beta; 280 281 do 282 { 283 cosTheta = 1. - 2. * G4UniformRand(); 284 fcos = (1 + cosTheta * cosTheta) * 285 } 286 // Loop checking, 07-Aug-2015, Vladimir 287 while(fcos < G4UniformRand()); 288 289 beta = std::sqrt(1. - 1. / (gamma * gamm 290 291 cosTheta = (cosTheta + beta) / (1. + bet 292 293 if(cosTheta > 1.) 294 cosTheta = 1.; 295 if(cosTheta < -1.) 296 cosTheta = -1.; 297 298 sinTheta = std::sqrt(1. - cosTheta * cos 299 300 G4double Phi = twopi * G4UniformRand(); 301 302 G4double dirx = sinTheta * std::cos(Phi) 303 G4double diry = sinTheta * std::sin(Phi) 304 G4double dirz = cosTheta; 305 306 G4ThreeVector gammaDirection(dirx, diry, 307 gammaDirection.rotateUz(particleDirectio 308 309 // polarization of new gamma 310 G4ThreeVector gammaPolarization = FieldV 311 gammaPolarization = gammaP 312 313 // create G4DynamicParticle object for t 314 auto aGamma = 315 new G4DynamicParticle(G4Gamma::Gamma() 316 aGamma->SetPolarization(gammaPolarizatio 317 gammaPolarizatio 318 319 aParticleChange.SetNumberOfSecondaries(1 320 321 // Update the incident particle 322 G4double newKinEnergy = kineticEnergy - 323 324 if(newKinEnergy > 0.) 325 { 326 aParticleChange.ProposeMomentumDirecti 327 aParticleChange.ProposeEnergy(newKinEn 328 aParticleChange.ProposeLocalEnergyDepo 329 } 330 else 331 { 332 aParticleChange.ProposeEnergy(0.); 333 aParticleChange.ProposeLocalEnergyDepo 334 G4double charge = aDynamicParticle->Ge 335 if(charge < 0.) 336 { 337 aParticleChange.ProposeTrackStatus(f 338 } 339 else 340 { 341 aParticleChange.ProposeTrackStatus(f 342 } 343 } 344 345 // Create the G4Track 346 G4Track* aSecondaryTrack = new G4Track(a 347 aSecondaryTrack->SetTouchableHandle(step 348 aSecondaryTrack->SetParentID(trackData.G 349 aSecondaryTrack->SetCreatorModelID(secID 350 aParticleChange.AddSecondary(aSecondaryT 351 } 352 else 353 { 354 return G4VDiscreteProcess::PostStepDoIt( 355 } 356 } 357 return G4VDiscreteProcess::PostStepDoIt(trac 358 } 359 360 G4double G4SynchrotronRadiationInMat::GetPhoto 361 362 { 363 G4int i; 364 G4double energyOfSR = -1.0; 365 366 const G4DynamicParticle* aDynamicParticle = 367 368 G4double gamma = 369 aDynamicParticle->GetTotalEnergy() / (aDyn 370 371 G4double particleCharge = aDynamicParticle-> 372 373 G4ThreeVector FieldValue; 374 const G4Field* pField = nullptr; 375 376 G4FieldManager* fieldMgr = nullptr; 377 G4bool fieldExertsForce = false; 378 379 if((particleCharge != 0.0)) 380 { 381 fieldMgr = fFieldPropagator->FindAndSetFie 382 if(fieldMgr != nullptr) 383 { 384 // If the field manager has no field, th 385 fieldExertsForce = (fieldMgr->GetDetecto 386 } 387 } 388 if(fieldExertsForce) 389 { 390 pField = fieldMgr->Get 391 G4ThreeVector globPosition = trackData.Get 392 G4double globPosVec[3], FieldValueVec[3]; 393 394 globPosVec[0] = globPosition.x(); 395 globPosVec[1] = globPosition.y(); 396 globPosVec[2] = globPosition.z(); 397 398 pField->GetFieldValue(globPosVec, FieldVal 399 FieldValue = 400 G4ThreeVector(FieldValueVec[0], FieldVal 401 402 G4ThreeVector unitMomentum = aDynamicParti 403 G4ThreeVector unitMcrossB = FieldValue.cr 404 G4double perpB = unitMcrossB.m 405 if(perpB > 0.0) 406 { 407 // M-C of synchrotron photon energy 408 G4double random = G4UniformRand(); 409 for(i = 0; i < 200; ++i) 410 { 411 if(random >= fIntegralProbabilityOfSR[ 412 break; 413 } 414 energyOfSR = 0.0001 * i * i * fEnergyCon 415 416 // check against insufficient energy 417 if(energyOfSR <= 0.0) 418 { 419 return -1.0; 420 } 421 } 422 else 423 { 424 return -1.0; 425 } 426 } 427 return energyOfSR; 428 } 429 430 ////////////////////////////////////////////// 431 G4double G4SynchrotronRadiationInMat::GetRando 432 433 { 434 G4int i; 435 static constexpr G4int iMax = 200; 436 G4double energySR, random, position; 437 438 random = G4UniformRand(); 439 440 for(i = 0; i < iMax; ++i) 441 { 442 if(random >= fIntegralProbabilityOfSR[i]) 443 break; 444 } 445 if(i <= 0) 446 position = G4UniformRand(); 447 else if(i >= iMax) 448 position = G4double(iMax); 449 else 450 position = i + G4UniformRand(); 451 452 energySR = 453 0.0001 * position * position * fEnergyCons 454 455 if(energySR < 0.) 456 energySR = 0.; 457 458 return energySR; 459 } 460 461 ////////////////////////////////////////////// 462 G4double G4SynchrotronRadiationInMat::GetProbS 463 { 464 G4double result, hypCos2, hypCos = std::cosh 465 466 hypCos2 = hypCos * hypCos; 467 result = std::cosh(5. * t / 3.) * std::exp(t 468 result /= hypCos2; 469 return result; 470 } 471 472 ////////////////////////////////////////////// 473 // return the probability to emit SR photon wi 474 // energy/energy_c >= ksi 475 // for ksi <= 0. P = 1., however the method wo 476 G4double G4SynchrotronRadiationInMat::GetIntPr 477 { 478 if(ksi <= 0.) 479 return 1.0; 480 fKsi = ksi; // should be > 0. ! 481 G4int n; 482 G4double result, a; 483 484 a = fAlpha; // always = 0. 485 n = fRootNumber; // around default = 80 486 487 G4Integrator<G4SynchrotronRadiationInMat, 488 G4double (G4SynchrotronRadiatio 489 integral; 490 491 result = integral.Laguerre( 492 this, &G4SynchrotronRadiationInMat::GetPro 493 494 result *= 3. / 5. / pi; 495 496 return result; 497 } 498 499 ////////////////////////////////////////////// 500 // return an auxiliary function for K_5/3 inte 501 G4double G4SynchrotronRadiationInMat::GetProbS 502 { 503 G4double result, hypCos = std::cosh(t); 504 505 result = std::cosh(5. * t / 3.) * std::exp(t 506 result /= hypCos; 507 return result; 508 } 509 510 ////////////////////////////////////////////// 511 // return the probability to emit SR photon en 512 // energy/energy_c >= ksi 513 // for ksi <= 0. P = 1., however the method wo 514 G4double G4SynchrotronRadiationInMat::GetEnerg 515 { 516 if(ksi <= 0.) 517 return 1.0; 518 fKsi = ksi; // should be > 0. ! 519 G4int n; 520 G4double result, a; 521 522 a = fAlpha; // always = 0. 523 n = fRootNumber; // around default = 80 524 525 G4Integrator<G4SynchrotronRadiationInMat, 526 G4double (G4SynchrotronRadiatio 527 integral; 528 529 result = integral.Laguerre( 530 this, &G4SynchrotronRadiationInMat::GetPro 531 532 result *= 9. * std::sqrt(3.) * ksi / 8. / pi 533 534 return result; 535 } 536 537 ////////////////////////////////////////////// 538 G4double G4SynchrotronRadiationInMat::GetInteg 539 { 540 G4double result, hypCos = std::cosh(t); 541 542 result = 543 std::cosh(fOrderAngleK * t) * std::exp(t - 544 result /= hypCos; 545 return result; 546 } 547 548 ////////////////////////////////////////////// 549 // Return K 1/3 or 2/3 for angular distributio 550 G4double G4SynchrotronRadiationInMat::GetAngle 551 { 552 fEta = eta; // should be > 0. ! 553 G4int n; 554 G4double result, a; 555 556 a = fAlpha; // always = 0. 557 n = fRootNumber; // around default = 80 558 559 G4Integrator<G4SynchrotronRadiationInMat, 560 G4double (G4SynchrotronRadiatio 561 integral; 562 563 result = integral.Laguerre( 564 this, &G4SynchrotronRadiationInMat::GetInt 565 566 return result; 567 } 568 569 ////////////////////////////////////////////// 570 // Relative angle diff distribution for given 571 G4double G4SynchrotronRadiationInMat::GetAngle 572 { 573 G4double result, funK, funK2, gpsi2 = gpsi * 574 575 fPsiGamma = gpsi; 576 fEta = 0.5 * fKsi * (1. + gpsi2) * std: 577 578 fOrderAngleK = 1. / 3.; 579 funK = GetAngleK(fEta); 580 funK2 = funK * funK; 581 582 result = gpsi2 * funK2 / (1. + gpsi2); 583 584 fOrderAngleK = 2. / 3.; 585 funK = GetAngleK(fEta); 586 funK2 = funK * funK; 587 588 result += funK2; 589 result *= (1. + gpsi2) * fKsi; 590 591 return result; 592 } 593