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 // Authors: G.Depaola & F.Longo 28 // 29 // History: 30 // ------- 31 // 32 // 05 Apr 2021 J Allison added quantum entan 33 // If the photons have been "tagged" as "quant 34 // G4eplusAnnihilation for annihilation into 2 35 // here if - and only if - both photons suffer 36 // predictions from Pryce and Ward, Nature No 37 // Physical Review 73 (1948) p.440. Experiment 38 // entanglement in the MeV regime and its appl 39 // D. Watts, J. Allison et al., Nature Communi 40 // https://doi.org/10.1038/s41467-021-22907-5. 41 // 42 // 02 May 2009 S Incerti as V. Ivanchenko pr 43 // 44 // Cleanup initialisation and generation of se 45 // - apply internal high-ener 46 // - do not apply low-energy 47 // - remove GetMeanFreePath m 48 // - added protection against 49 // - use G4ElementSelector 50 51 #include "G4LivermorePolarizedComptonModel.hh" 52 #include "G4PhysicalConstants.hh" 53 #include "G4SystemOfUnits.hh" 54 #include "G4AutoLock.hh" 55 #include "G4Electron.hh" 56 #include "G4ParticleChangeForGamma.hh" 57 #include "G4LossTableManager.hh" 58 #include "G4VAtomDeexcitation.hh" 59 #include "G4AtomicShell.hh" 60 #include "G4Gamma.hh" 61 #include "G4ShellData.hh" 62 #include "G4DopplerProfile.hh" 63 #include "G4Log.hh" 64 #include "G4Exp.hh" 65 #include "G4Pow.hh" 66 #include "G4LogLogInterpolation.hh" 67 #include "G4PhysicsModelCatalog.hh" 68 #include "G4EntanglementAuxInfo.hh" 69 #include "G4eplusAnnihilationEntanglementClipB 70 71 //....oooOO0OOooo........oooOO0OOooo........oo 72 73 using namespace std; 74 namespace { G4Mutex LivermorePolarizedComptonM 75 76 77 G4PhysicsFreeVector* G4LivermorePolarizedCompt 78 G4ShellData* G4LivermorePolarizedCompton 79 G4DopplerProfile* G4LivermorePolarizedCompton 80 G4CompositeEMDataSet* G4LivermorePolarizedComp 81 82 //....oooOO0OOooo........oooOO0OOooo........oo 83 84 G4LivermorePolarizedComptonModel::G4LivermoreP 85 :G4VEmModel(nam),isInitialised(false) 86 { 87 verboseLevel= 1; 88 // Verbosity scale: 89 // 0 = nothing 90 // 1 = warning for energy non-conservation 91 // 2 = details of energy budget 92 // 3 = calculation of cross sections, file o 93 // 4 = entering in methods 94 95 if( verboseLevel>1 ) 96 G4cout << "Livermore Polarized Compton is 97 98 //Mark this model as "applicable" for atomic 99 SetDeexcitationFlag(true); 100 101 fParticleChange = nullptr; 102 fAtomDeexcitation = nullptr; 103 fEntanglementModelID = G4PhysicsModelCatalog 104 } 105 106 //....oooOO0OOooo........oooOO0OOooo........oo 107 108 G4LivermorePolarizedComptonModel::~G4Livermore 109 { 110 if(IsMaster()) { 111 delete shellData; 112 shellData = nullptr; 113 delete profileData; 114 profileData = nullptr; 115 delete scatterFunctionData; 116 scatterFunctionData = nullptr; 117 for(G4int i=0; i<maxZ; ++i) { 118 if(data[i]) { 119 delete data[i]; 120 data[i] = nullptr; 121 } 122 } 123 } 124 } 125 126 //....oooOO0OOooo........oooOO0OOooo........oo 127 128 void G4LivermorePolarizedComptonModel::Initial 129 const G 130 { 131 if (verboseLevel > 1) 132 G4cout << "Calling G4LivermorePolarizedCom 133 134 // Initialise element selector 135 if(IsMaster()) { 136 // Access to elements 137 const char* path = G4FindDataDir("G4LEDATA 138 139 G4ProductionCutsTable* theCoupleTable = 140 G4ProductionCutsTable::GetProductionCuts 141 142 G4int numOfCouples = (G4int)theCoupleTable 143 144 for(G4int i=0; i<numOfCouples; ++i) { 145 const G4Material* material = 146 theCoupleTable->GetMaterialCutsCouple( 147 const G4ElementVector* theElementVector 148 std::size_t nelm = material->GetNumberOf 149 150 for (std::size_t j=0; j<nelm; ++j) { 151 G4int Z = G4lrint((*theElementVector)[ 152 if(Z < 1) { Z = 1; } 153 else if(Z > maxZ){ Z = maxZ; } 154 155 if( (!data[Z]) ) { ReadData(Z, path); 156 } 157 } 158 159 // For Doppler broadening 160 if(!shellData) { 161 shellData = new G4ShellData(); 162 shellData->SetOccupancyData(); 163 G4String file = "/doppler/shell-doppler" 164 shellData->LoadData(file); 165 } 166 if(!profileData) { profileData = new G4Dop 167 168 // Scattering Function 169 if(!scatterFunctionData) 170 { 171 172 G4VDataSetAlgorithm* scatterInterpolation = 173 G4String scatterFile = "comp/ce-sf-"; 174 scatterFunctionData = new G4CompositeEMDataS 175 scatterFunctionData->LoadData(scatterFile); 176 } 177 178 InitialiseElementSelectors(particle, cuts) 179 } 180 181 if (verboseLevel > 2) { 182 G4cout << "Loaded cross section files" << 183 } 184 185 if( verboseLevel>1 ) { 186 G4cout << "G4LivermoreComptonModel is init 187 << "Energy range: " 188 << LowEnergyLimit() / eV << " eV - " 189 << HighEnergyLimit() / GeV << " GeV" 190 << G4endl; 191 } 192 // 193 if(isInitialised) { return; } 194 195 fParticleChange = GetParticleChangeForGamma( 196 fAtomDeexcitation = G4LossTableManager::Ins 197 isInitialised = true; 198 } 199 200 201 void G4LivermorePolarizedComptonModel::Initial 202 G4VEmModel* masterModel) 203 { 204 SetElementSelectors(masterModel->GetElementS 205 } 206 207 //....oooOO0OOooo........oooOO0OOooo........oo 208 209 void G4LivermorePolarizedComptonModel::ReadDat 210 { 211 if (verboseLevel > 1) 212 { 213 G4cout << "G4LivermorePolarizedComptonMo 214 << G4endl; 215 } 216 if(data[Z]) { return; } 217 const char* datadir = path; 218 if(!datadir) 219 { 220 datadir = G4FindDataDir("G4LEDATA"); 221 if(!datadir) 222 { 223 G4Exception("G4LivermorePolarizedComptonMo 224 "em0006",FatalException, 225 "Environment variable G4LEDATA not d 226 return; 227 } 228 } 229 230 data[Z] = new G4PhysicsFreeVector(); 231 232 std::ostringstream ost; 233 ost << datadir << "/livermore/comp/ce-cs-" < 234 std::ifstream fin(ost.str().c_str()); 235 236 if( !fin.is_open()) 237 { 238 G4ExceptionDescription ed; 239 ed << "G4LivermorePolarizedComptonModel 240 << "> is not opened!" << G4endl; 241 G4Exception("G4LivermoreComptonModel::Re 242 "em0003",FatalException, 243 ed,"G4LEDATA version should be G4EMLOW8. 244 return; 245 } else { 246 if(verboseLevel > 3) { 247 G4cout << "File " << ost.str() 248 << " is opened by G4LivermorePolarizedC 249 } 250 data[Z]->Retrieve(fin, true); 251 data[Z]->ScaleVector(MeV, MeV*barn); 252 } 253 fin.close(); 254 } 255 256 //....oooOO0OOooo........oooOO0OOooo........oo 257 258 G4double G4LivermorePolarizedComptonModel::Com 259 const G 260 G 261 G 262 G 263 { 264 if (verboseLevel > 3) 265 G4cout << "Calling ComputeCrossSectionPerA 266 267 G4double cs = 0.0; 268 269 if (GammaEnergy < LowEnergyLimit()) 270 return 0.0; 271 272 G4int intZ = G4lrint(Z); 273 if(intZ < 1 || intZ > maxZ) { return cs; } 274 275 G4PhysicsFreeVector* pv = data[intZ]; 276 277 // if element was not initialised 278 // do initialisation safely for MT mode 279 if(!pv) 280 { 281 InitialiseForElement(0, intZ); 282 pv = data[intZ]; 283 if(!pv) { return cs; } 284 } 285 286 G4int n = G4int(pv->GetVectorLength() - 1); 287 G4double e1 = pv->Energy(0); 288 G4double e2 = pv->Energy(n); 289 290 if(GammaEnergy <= e1) { cs = GammaEnerg 291 else if(GammaEnergy <= e2) { cs = pv->Value( 292 else if(GammaEnergy > e2) { cs = pv->Value( 293 294 return cs; 295 296 } 297 298 //....oooOO0OOooo........oooOO0OOooo........oo 299 300 void G4LivermorePolarizedComptonModel::SampleS 301 const G4MaterialCutsCouple* co 302 const G4DynamicParticle* aDyna 303 G4double, 304 G4double) 305 { 306 // The scattered gamma energy is sampled acc 307 // The random number techniques of Butcher & 308 // GEANT4 internal units 309 // 310 // Note : Effects due to binding of atomic e 311 312 if (verboseLevel > 3) 313 G4cout << "Calling SampleSecondaries() of 314 315 G4double gammaEnergy0 = aDynamicGamma->GetKi 316 317 // do nothing below the threshold 318 // should never get here because the XS is z 319 if (gammaEnergy0 < LowEnergyLimit()) 320 return ; 321 322 G4ThreeVector gammaPolarization0 = aDynamicG 323 324 // Protection: a polarisation parallel to th 325 // direction causes problems; 326 // in that case find a random polarization 327 G4ThreeVector gammaDirection0 = aDynamicGamm 328 329 // Make sure that the polarization vector is 330 // gamma direction. If not 331 if(!(gammaPolarization0.isOrthogonal(gammaDi 332 { // only for testing now 333 gammaPolarization0 = GetRandomPolarizati 334 } 335 else 336 { 337 if ( gammaPolarization0.howOrthogonal(ga 338 { 339 gammaPolarization0 = GetPerpendicularPolar 340 } 341 } 342 // End of Protection 343 344 G4double E0_m = gammaEnergy0 / electron_mass 345 346 // Select randomly one element in the curren 347 //G4int Z = crossSectionHandler->SelectRando 348 const G4ParticleDefinition* particle = aDyn 349 const G4Element* elm = SelectRandomAtom(coup 350 G4int Z = (G4int)elm->GetZ(); 351 352 // Sample the energy and the polarization of 353 G4double epsilon, epsilonSq, onecost, sinThe 354 355 G4double epsilon0Local = 1./(1. + 2*E0_m); 356 G4double epsilon0Sq = epsilon0Local*epsilon0 357 G4double alpha1 = - G4Log(epsilon0Local); 358 G4double alpha2 = 0.5*(1.- epsilon0Sq); 359 360 G4double wlGamma = h_Planck*c_light/gammaEne 361 G4double gammaEnergy1; 362 G4ThreeVector gammaDirection1; 363 364 do { 365 if ( alpha1/(alpha1+alpha2) > G4UniformRan 366 { 367 epsilon = G4Exp(-alpha1*G4UniformRand()); 368 epsilonSq = epsilon*epsilon; 369 } 370 else 371 { 372 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4 373 epsilon = std::sqrt(epsilonSq); 374 } 375 376 onecost = (1.- epsilon)/(epsilon*E0_m); 377 sinThetaSqr = onecost*(2.-onecost); 378 379 // Protection 380 if (sinThetaSqr > 1.) 381 { 382 G4cout 383 << " -- Warning -- G4LivermorePolarizedCom 384 << "sin(theta)**2 = " 385 << sinThetaSqr 386 << "; set to 1" 387 << G4endl; 388 sinThetaSqr = 1.; 389 } 390 if (sinThetaSqr < 0.) 391 { 392 G4cout 393 << " -- Warning -- G4LivermorePolarizedCom 394 << "sin(theta)**2 = " 395 << sinThetaSqr 396 << "; set to 0" 397 << G4endl; 398 sinThetaSqr = 0.; 399 } 400 // End protection 401 402 G4double x = std::sqrt(onecost/2.) / (wlG 403 G4double scatteringFunction = scatterFunct 404 greject = (1. - epsilon*sinThetaSqr/(1.+ e 405 406 } while(greject < G4UniformRand()*Z); 407 408 409 // ***************************************** 410 // Phi determination 411 // ***************************************** 412 G4double phi = SetPhi(epsilon,sinThetaSqr); 413 414 // 415 // scattered gamma angles. ( Z - axis along 416 // 417 G4double cosTheta = 1. - onecost; 418 419 // Protection 420 if (cosTheta > 1.) 421 { 422 G4cout 423 << " -- Warning -- G4LivermorePolarizedCompt 424 << "cosTheta = " 425 << cosTheta 426 << "; set to 1" 427 << G4endl; 428 cosTheta = 1.; 429 } 430 if (cosTheta < -1.) 431 { 432 G4cout 433 << " -- Warning -- G4LivermorePolarizedCompt 434 << "cosTheta = " 435 << cosTheta 436 << "; set to -1" 437 << G4endl; 438 cosTheta = -1.; 439 } 440 // End protection 441 442 G4double sinTheta = std::sqrt (sinThetaSqr); 443 444 // Protection 445 if (sinTheta > 1.) 446 { 447 G4cout 448 << " -- Warning -- G4LivermorePolarizedCompt 449 << "sinTheta = " 450 << sinTheta 451 << "; set to 1" 452 << G4endl; 453 sinTheta = 1.; 454 } 455 if (sinTheta < -1.) 456 { 457 G4cout 458 << " -- Warning -- G4LivermorePolarizedCompt 459 << "sinTheta = " 460 << sinTheta 461 << "; set to -1" 462 << G4endl; 463 sinTheta = -1.; 464 } 465 // End protection 466 467 // Check for entanglement and re-sample phi 468 469 const auto* auxInfo 470 = fParticleChange->GetCurrentTrack()->GetAux 471 if (auxInfo) { 472 const auto* entanglementAuxInfo = dynamic_ 473 if (entanglementAuxInfo) { 474 auto* clipBoard = dynamic_cast<G4eplusAn 475 (entanglementAuxInfo->GetEntanglementCli 476 if (clipBoard) { 477 // This is an entangled photon from ep 478 // If this is the first scatter of the 479 // phi on the clipboard. 480 // If this is the first scatter of the 481 // phi of the first scatter of the fir 482 // theta of the second photon, to samp 483 if (clipBoard->IsTrack1Measurement()) 484 // Check we have the relevant track. 485 // necessary but I want to be sure t 486 // entangled system are properly pai 487 // Note: the tracking manager pops t 488 // will rely on that. (If not, the 489 // more complicated to ensure we mat 490 // So our track 1 is clipboard track 491 if (clipBoard->GetTrackB() == fParti 492 // This is the first scatter of th 493 // // Debug 494 // auto* track1 = fPart 495 // G4cout 496 // << "This is the firs 497 // << "\nTrack: " << tr 498 // << ", Parent: " << t 499 // << ", Name: " << cli 500 // << G4endl; 501 // // End debug 502 clipBoard->ResetTrack1Measurement( 503 // Store cos(theta),phi of first p 504 clipBoard->SetComptonCosTheta1(cos 505 clipBoard->SetComptonPhi1(phi); 506 } 507 } else if (clipBoard->IsTrack2Measurem 508 // Check we have the relevant track. 509 // Remember our track 2 is clipboard 510 if (clipBoard->GetTrackA() == fParti 511 // This is the first scatter of th 512 // // Debug 513 // auto* track2 = fPart 514 // G4cout 515 // << "This is the firs 516 // << "\nTrack: " << tr 517 // << ", Parent: " << t 518 // << ", Name: " << cli 519 // << G4endl; 520 // // End debug 521 clipBoard->ResetTrack2Measurement( 522 523 // Get cos(theta),phi of first pho 524 const G4double& cosTheta1 = clipBo 525 const G4double& phi1 = clipBoard-> 526 // For clarity make aliases for th 527 const G4double& cosTheta2 = cosThe 528 G4double& phi2 = phi; 529 // G4cout << "cosTheta1 530 // G4cout << "cosTheta2 531 532 // Re-sample phi 533 // Draw the difference of azimutha 534 // A + B * cos(2*deltaPhi), or rat 535 // C = A / (A + |B|) and D = B / ( 536 const G4double sin2Theta1 = 1.-cos 537 const G4double sin2Theta2 = 1.-cos 538 539 // Pryce and Ward, Nature No 4065 540 auto* g4Pow = G4Pow::GetInstance() 541 const G4double A = 542 ((g4Pow->powN(1.-cosTheta1,3))+2.) 543 ((g4Pow->powN(2.-cosTheta1,3)*g4Po 544 const G4double B = -(sin2Theta1*si 545 ((g4Pow->powN(2.-cosTheta1,2)*g4Po 546 547 // // Snyder et al, Physical Re 548 // // (This is an alternative f 549 // const G4double& k0 = gammaEn 550 // const G4double k1 = k0/(2.-c 551 // const G4double k2 = k0/(2.-c 552 // const G4double gamma1 = k1/k 553 // const G4double gamma2 = k2/k 554 // const G4double A1 = gamma1*g 555 // const G4double B1 = 2.*sin2T 556 // // That's A1 + B1*sin2(delta 557 // const G4double A = A1 + 0.5* 558 // const G4double B = -0.5*B1; 559 560 const G4double maxValue = A + std: 561 const G4double C = A / maxValue; 562 const G4double D = B / maxValue; 563 // G4cout << "A,B,C,D: " << A < 564 565 // Sample delta phi 566 G4double deltaPhi; 567 const G4int maxCount = 999999; 568 G4int iCount = 0; 569 for (; iCount < maxCount; ++iCount 570 deltaPhi = twopi * G4UniformRand 571 if (G4UniformRand() < C + D * co 572 } 573 if (iCount >= maxCount ) { 574 G4cout << "G4LivermorePolarizedC 575 << "Re-sampled delta phi not fou 576 << " tries - carrying on anyway. 577 } 578 579 // Thus, the desired second photon 580 phi2 = deltaPhi - phi1 + halfpi; 581 // The minus sign is in above stat 582 // annihilation photons are in opp 583 // are measured in the opposite di 584 // halfpi is added for the followi 585 // In this function phi is relativ 586 // SystemOfRefChange below. We kno 587 // the polarisations of the two an 588 // to each other, i.e., halfpi dif 589 // Furthermore, only sin(phi) and 590 // need to place any range constra 591 // if (phi2 > pi) { 592 // phi2 -= twopi; 593 // } 594 // if (phi2 < -pi) { 595 // phi2 += twopi; 596 // } 597 } 598 } 599 } 600 } 601 } 602 603 // End of entanglement 604 605 G4double dirx = sinTheta*std::cos(phi); 606 G4double diry = sinTheta*std::sin(phi); 607 G4double dirz = cosTheta ; 608 609 // oneCosT , eom 610 611 // Doppler broadening - Method based on: 612 // Y. Namito, S. Ban and H. Hirayama, 613 // "Implementation of the Doppler Broadening 614 // NIM A 349, pp. 489-494, 1994 615 616 // Maximum number of sampling iterations 617 static G4int maxDopplerIterations = 1000; 618 G4double bindingE = 0.; 619 G4double photonEoriginal = epsilon * gammaEn 620 G4double photonE = -1.; 621 G4int iteration = 0; 622 G4double eMax = gammaEnergy0; 623 624 G4int shellIdx = 0; 625 626 if (verboseLevel > 3) { 627 G4cout << "Started loop to sample broading 628 } 629 630 do 631 { 632 iteration++; 633 // Select shell based on shell occupancy 634 shellIdx = shellData->SelectRandomShell( 635 bindingE = shellData->BindingEnergy(Z,sh 636 637 if (verboseLevel > 3) { 638 G4cout << "Shell ID= " << shellIdx 639 << " Ebind(keV)= " << bindingE/keV << 640 } 641 eMax = gammaEnergy0 - bindingE; 642 643 // Randomly sample bound electron moment 644 G4double pSample = profileData->RandomSe 645 646 if (verboseLevel > 3) { 647 G4cout << "pSample= " << pSample << G4e 648 } 649 // Rescale from atomic units 650 G4double pDoppler = pSample * fine_struc 651 G4double pDoppler2 = pDoppler * pDoppler 652 G4double var2 = 1. + onecost * E0_m; 653 G4double var3 = var2*var2 - pDoppler2; 654 G4double var4 = var2 - pDoppler2 * cosTh 655 G4double var = var4*var4 - var3 + pDoppl 656 if (var > 0.) 657 { 658 G4double varSqrt = std::sqrt(var); 659 G4double scale = gammaEnergy0 / var3; 660 // Random select either root 661 if (G4UniformRand() < 0.5) photonE = (var4 662 else photonE = (var4 + varSqrt) * scale; 663 } 664 else 665 { 666 photonE = -1.; 667 } 668 } while ( iteration <= maxDopplerIterations 669 (photonE < 0. || photonE > eMax || phot 670 671 // End of recalculation of photon energy wit 672 // Revert to original if maximum number of i 673 if (iteration >= maxDopplerIterations) 674 { 675 photonE = photonEoriginal; 676 bindingE = 0.; 677 } 678 679 gammaEnergy1 = photonE; 680 681 // 682 // update G4VParticleChange for the scattere 683 // 684 // New polarization 685 G4ThreeVector gammaPolarization1 = SetNewPol 686 sinThetaSqr, 687 phi, 688 cosTheta); 689 690 // Set new direction 691 G4ThreeVector tmpDirection1( dirx,diry,dirz 692 gammaDirection1 = tmpDirection1; 693 694 // Change reference frame. 695 SystemOfRefChange(gammaDirection0,gammaDirec 696 gammaPolarization0,gammaPolarization1) 697 698 if (gammaEnergy1 > 0.) 699 { 700 fParticleChange->SetProposedKineticEnerg 701 fParticleChange->ProposeMomentumDirectio 702 fParticleChange->ProposePolarization( ga 703 } 704 else 705 { 706 gammaEnergy1 = 0.; 707 fParticleChange->SetProposedKineticEnerg 708 fParticleChange->ProposeTrackStatus(fSto 709 } 710 711 // 712 // kinematic of the scattered electron 713 // 714 G4double ElecKineEnergy = gammaEnergy0 - gam 715 716 // SI -protection against negative final ene 717 // like in G4LivermoreComptonModel.cc 718 if(ElecKineEnergy < 0.0) { 719 fParticleChange->ProposeLocalEnergyDeposit 720 return; 721 } 722 723 G4double ElecMomentum = std::sqrt(ElecKineEn 724 725 G4ThreeVector ElecDirection((gammaEnergy0 * 726 gammaEnergy1 * gammaDirection1) * ( 727 728 G4DynamicParticle* dp = 729 new G4DynamicParticle (G4Electron::Electro 730 fvect->push_back(dp); 731 732 // sample deexcitation 733 // 734 if (verboseLevel > 3) { 735 G4cout << "Started atomic de-excitation " 736 } 737 738 if(fAtomDeexcitation && iteration < maxDoppl 739 G4int index = couple->GetIndex(); 740 if(fAtomDeexcitation->CheckDeexcitationAct 741 std::size_t nbefore = fvect->size(); 742 G4AtomicShellEnumerator as = G4AtomicShe 743 const G4AtomicShell* shell = fAtomDeexci 744 fAtomDeexcitation->GenerateParticles(fve 745 std::size_t nafter = fvect->size(); 746 if(nafter > nbefore) { 747 for (std::size_t i=nbefore; i<nafter; ++i) { 748 //Check if there is enough residual energy 749 if (bindingE >= ((*fvect)[i])->GetKineticE 750 { 751 //Ok, this is a valid secondary: 752 bindingE -= ((*fvect)[i])->GetKineticE 753 } 754 else 755 { 756 //Invalid secondary: not enough energy 757 //Keep its energy in the local deposit 758 delete (*fvect)[i]; 759 (*fvect)[i]=0; 760 } 761 } 762 } 763 } 764 } 765 //This should never happen 766 if(bindingE < 0.0) 767 G4Exception("G4LivermoreComptonModel::Samp 768 "em2050",FatalException,"Negative local en 769 770 fParticleChange->ProposeLocalEnergyDeposit(b 771 772 } 773 774 //....oooOO0OOooo........oooOO0OOooo........oo 775 776 G4double G4LivermorePolarizedComptonModel::Set 777 G4double sinSqrTh) 778 { 779 G4double rand1; 780 G4double rand2; 781 G4double phiProbability; 782 G4double phi; 783 G4double a, b; 784 785 do 786 { 787 rand1 = G4UniformRand(); 788 rand2 = G4UniformRand(); 789 phiProbability=0.; 790 phi = twopi*rand1; 791 792 a = 2*sinSqrTh; 793 b = energyRate + 1/energyRate; 794 795 phiProbability = 1 - (a/b)*(std::cos(phi 796 } 797 while ( rand2 > phiProbability ); 798 return phi; 799 } 800 801 802 //....oooOO0OOooo........oooOO0OOooo........oo 803 804 G4ThreeVector G4LivermorePolarizedComptonModel 805 { 806 G4double dx = a.x(); 807 G4double dy = a.y(); 808 G4double dz = a.z(); 809 G4double x = dx < 0.0 ? -dx : dx; 810 G4double y = dy < 0.0 ? -dy : dy; 811 G4double z = dz < 0.0 ? -dz : dz; 812 if (x < y) { 813 return x < z ? G4ThreeVector(-dy,dx,0) : G 814 }else{ 815 return y < z ? G4ThreeVector(dz,0,-dx) : G 816 } 817 } 818 819 //....oooOO0OOooo........oooOO0OOooo........oo 820 821 G4ThreeVector G4LivermorePolarizedComptonModel 822 { 823 G4ThreeVector d0 = direction0.unit(); 824 G4ThreeVector a1 = SetPerpendicularVector(d0 825 G4ThreeVector a0 = a1.unit(); // unit vector 826 827 G4double rand1 = G4UniformRand(); 828 829 G4double angle = twopi*rand1; // random pola 830 G4ThreeVector b0 = d0.cross(a0); // cross pr 831 832 G4ThreeVector c; 833 834 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 835 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 836 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 837 838 G4ThreeVector c0 = c.unit(); 839 840 return c0; 841 } 842 843 //....oooOO0OOooo........oooOO0OOooo........oo 844 845 G4ThreeVector G4LivermorePolarizedComptonModel 846 (const G4ThreeVector& gammaDirection, const G4 847 { 848 // 849 // The polarization of a photon is always pe 850 // Therefore this function removes those vec 851 // points in direction of gammaDirection 852 // 853 // Mathematically we search the projection o 854 // plains normal vector. 855 // The basic equation can be found in each g 856 // p = a - (a o n)/(n o n)*n 857 858 return gammaPolarization - gammaPolarization 859 } 860 861 //....oooOO0OOooo........oooOO0OOooo........oo 862 863 G4ThreeVector G4LivermorePolarizedComptonModel 864 G4double sinSqrTh, 865 G4double phi, 866 G4double costheta) 867 { 868 G4double rand1; 869 G4double rand2; 870 G4double cosPhi = std::cos(phi); 871 G4double sinPhi = std::sin(phi); 872 G4double sinTheta = std::sqrt(sinSqrTh); 873 G4double cosSqrPhi = cosPhi*cosPhi; 874 // G4double cossqrth = 1.-sinSqrTh; 875 // G4double sinsqrphi = sinPhi*sinPhi; 876 G4double normalisation = std::sqrt(1. - cosS 877 878 // Determination of Theta 879 G4double theta; 880 881 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 882 rand1 = G4UniformRand(); 883 rand2 = G4UniformRand(); 884 885 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi 886 { 887 if (rand2<0.5) 888 theta = pi/2.0; 889 else 890 theta = 3.0*pi/2.0; 891 } 892 else 893 { 894 if (rand2<0.5) 895 theta = 0; 896 else 897 theta = pi; 898 } 899 G4double cosBeta = std::cos(theta); 900 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 901 902 G4ThreeVector gammaPolarization1; 903 904 G4double xParallel = normalisation*cosBeta; 905 G4double yParallel = -(sinSqrTh*cosPhi*sinPh 906 G4double zParallel = -(costheta*sinTheta*cos 907 G4double xPerpendicular = 0.; 908 G4double yPerpendicular = (costheta)*sinBeta 909 G4double zPerpendicular = -(sinTheta*sinPhi) 910 911 G4double xTotal = (xParallel + xPerpendicula 912 G4double yTotal = (yParallel + yPerpendicula 913 G4double zTotal = (zParallel + zPerpendicula 914 915 gammaPolarization1.setX(xTotal); 916 gammaPolarization1.setY(yTotal); 917 gammaPolarization1.setZ(zTotal); 918 919 return gammaPolarization1; 920 } 921 922 //....oooOO0OOooo........oooOO0OOooo........oo 923 924 void G4LivermorePolarizedComptonModel::SystemO 925 G4ThreeVector& direction1, 926 G4ThreeVector& polarization0, 927 G4ThreeVector& polarization1) 928 { 929 // direction0 is the original photon directi 930 // polarization0 is the original photon pola 931 // need to specify y axis in the real refere 932 G4ThreeVector Axis_Z0 = direction0.unit(); 933 G4ThreeVector Axis_X0 = polarization0.unit() 934 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 935 936 G4double direction_x = direction1.getX(); 937 G4double direction_y = direction1.getY(); 938 G4double direction_z = direction1.getZ(); 939 940 direction1 = (direction_x*Axis_X0 + directio 941 G4double polarization_x = polarization1.getX 942 G4double polarization_y = polarization1.getY 943 G4double polarization_z = polarization1.getZ 944 945 polarization1 = (polarization_x*Axis_X0 + po 946 947 } 948 949 950 //....oooOO0OOooo........oooOO0OOooo........oo 951 //....oooOO0OOooo........oooOO0OOooo........oo 952 953 void 954 G4LivermorePolarizedComptonModel::InitialiseFo 955 G4int Z) 956 { 957 G4AutoLock l(&LivermorePolarizedComptonModel 958 if(!data[Z]) { ReadData(Z); } 959 l.unlock(); 960 } 961