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 // | G4LowEPPolarizedComptonModel-- Geant 29 // | polarised low energy Compton scat 30 // | J. M. C. Brown, Monash Universit 31 // | 32 // | 33 // ******************************************* 34 // | 35 // | The following is a Geant4 class to simula 36 // | bound electron Compton scattering. Genera 37 // | based on G4LowEnergyCompton.cc and 38 // | G4LivermorePolarizedComptonModel.cc. 39 // | Algorithms for photon energy, and ejected 40 // | direction taken from: 41 // | 42 // | J. M. C. Brown, M. R. Dimmock, J. E. Gill 43 // | "A low energy bound atomic electron Compt 44 // | for Geant4", NIMB, Vol. 338, 77-88, 2014 45 // | 46 // | The author acknowledges the work of the G 47 // | in developing the following algorithms th 48 // | or adapeted for the present software: 49 // | 50 // | # sampling of photon scattering angle, 51 // | # target element selection in composite 52 // | # target shell selection in element, 53 // | # and sampling of bound electron momentu 54 // | 55 // ******************************************* 56 // | 57 // | History: 58 // | -------- 59 // | 60 // | Jan. 2015 JMCB - 1st Version based 61 // | Feb. 2016 JMCB - Geant4 10.2 FPE fi 62 // | Nov. 2016 JMCB - Polarisation track 63 // | of Dr. Merlin Reyn 64 // | University of Gene 65 // | 66 // ******************************************* 67 68 #include "G4LowEPPolarizedComptonModel.hh" 69 #include "G4AutoLock.hh" 70 #include "G4PhysicalConstants.hh" 71 #include "G4SystemOfUnits.hh" 72 73 //******************************************** 74 75 using namespace std; 76 namespace { G4Mutex LowEPPolarizedComptonModel 77 78 79 G4PhysicsFreeVector* G4LowEPPolarizedComptonMo 80 G4ShellData* G4LowEPPolarizedComptonMode 81 G4DopplerProfile* G4LowEPPolarizedComptonMode 82 83 static const G4double ln10 = G4Log(10.); 84 85 G4LowEPPolarizedComptonModel::G4LowEPPolarized 86 87 : G4VEmModel(nam),isInitialised(false) 88 { 89 verboseLevel=1 ; 90 // Verbosity scale: 91 // 0 = nothing 92 // 1 = warning for energy non-conservation 93 // 2 = details of energy budget 94 // 3 = calculation of cross sections, file o 95 // 4 = entering in methods 96 97 if( verboseLevel>1 ) { 98 G4cout << "Low energy photon Compton model 99 } 100 101 //Mark this model as "applicable" for atomic 102 SetDeexcitationFlag(true); 103 104 fParticleChange = nullptr; 105 fAtomDeexcitation = nullptr; 106 } 107 108 //******************************************** 109 110 G4LowEPPolarizedComptonModel::~G4LowEPPolarize 111 { 112 if(IsMaster()) { 113 delete shellData; 114 shellData = nullptr; 115 delete profileData; 116 profileData = nullptr; 117 } 118 } 119 120 //******************************************** 121 122 void G4LowEPPolarizedComptonModel::Initialise( 123 const 124 { 125 if (verboseLevel > 1) { 126 G4cout << "Calling G4LowEPPolarizedCompton 127 } 128 129 // Initialise element selector 130 if(IsMaster()) { 131 // Access to elements 132 const char* path = G4FindDataDir("G4LEDATA 133 134 G4ProductionCutsTable* theCoupleTable = 135 G4ProductionCutsTable::GetProductionCuts 136 G4int numOfCouples = (G4int)theCoupleTable 137 138 for(G4int i=0; i<numOfCouples; ++i) { 139 const G4Material* material = 140 theCoupleTable->GetMaterialCutsCouple( 141 const G4ElementVector* theElementVector 142 std::size_t nelm = material->GetNumberOf 143 144 for (std::size_t j=0; j<nelm; ++j) { 145 G4int Z = G4lrint((*theElementVector)[ 146 if(Z < 1) { Z = 1; } 147 else if(Z > maxZ){ Z = maxZ; } 148 149 if( (!data[Z]) ) { ReadData(Z, path); 150 } 151 } 152 153 // For Doppler broadening 154 if(!shellData) { 155 shellData = new G4ShellData(); 156 shellData->SetOccupancyData(); 157 G4String file = "/doppler/shell-doppler" 158 shellData->LoadData(file); 159 } 160 if(!profileData) { profileData = new G4Dop 161 162 InitialiseElementSelectors(particle, cuts) 163 } 164 165 if (verboseLevel > 2) { 166 G4cout << "Loaded cross section files" << 167 } 168 169 if( verboseLevel>1 ) { 170 G4cout << "G4LowEPPolarizedComptonModel is 171 << "Energy range: " 172 << LowEnergyLimit() / eV << " eV - 173 << HighEnergyLimit() / GeV << " GeV 174 << G4endl; 175 } 176 177 if(isInitialised) { return; } 178 179 fParticleChange = GetParticleChangeForGamma( 180 fAtomDeexcitation = G4LossTableManager::Ins 181 isInitialised = true; 182 } 183 184 //******************************************** 185 186 void G4LowEPPolarizedComptonModel::InitialiseL 187 188 { 189 SetElementSelectors(masterModel->GetElementS 190 } 191 192 //******************************************** 193 194 void G4LowEPPolarizedComptonModel::ReadData(st 195 { 196 if (verboseLevel > 1) 197 { 198 G4cout << "G4LowEPPolarizedComptonModel::R 199 << G4endl; 200 } 201 if(data[Z]) { return; } 202 const char* datadir = path; 203 if(!datadir) 204 { 205 datadir = G4FindDataDir("G4LEDATA"); 206 if(!datadir) 207 { 208 G4Exception("G4LowEPPolarizedComptonMode 209 "em0006",FatalException, 210 "Environment variable G4LEDA 211 return; 212 } 213 } 214 215 data[Z] = new G4PhysicsFreeVector(); 216 217 std::ostringstream ost; 218 ost << datadir << "/livermore/comp/ce-cs-" < 219 std::ifstream fin(ost.str().c_str()); 220 221 if( !fin.is_open()) 222 { 223 G4ExceptionDescription ed; 224 ed << "G4LowEPPolarizedComptonModel data 225 << "> is not opened!" << G4endl; 226 G4Exception("G4LowEPPolarizedComptonModel: 227 "em0003",FatalException, 228 ed,"G4LEDATA version should be 229 return; 230 } else { 231 if(verboseLevel > 3) { 232 G4cout << "File " << ost.str() 233 << " is opened by G4LowEPPolari 234 } 235 data[Z]->Retrieve(fin, true); 236 data[Z]->ScaleVector(MeV, MeV*barn); 237 } 238 fin.close(); 239 } 240 241 //******************************************** 242 243 G4double 244 G4LowEPPolarizedComptonModel::ComputeCrossSect 245 246 247 248 { 249 if (verboseLevel > 3) { 250 G4cout << "G4LowEPPolarizedComptonModel::C 251 << G4endl; 252 } 253 G4double cs = 0.0; 254 255 if (GammaEnergy < LowEnergyLimit()) { return 256 257 G4int intZ = G4lrint(Z); 258 if(intZ < 1 || intZ > maxZ) { return cs; } 259 260 G4PhysicsFreeVector* pv = data[intZ]; 261 262 // if element was not initialised 263 // do initialisation safely for MT mode 264 if(!pv) 265 { 266 InitialiseForElement(0, intZ); 267 pv = data[intZ]; 268 if(!pv) { return cs; } 269 } 270 271 G4int n = G4int(pv->GetVectorLength() - 1); 272 G4double e1 = pv->Energy(0); 273 G4double e2 = pv->Energy(n); 274 275 if(GammaEnergy <= e1) { cs = GammaEnerg 276 else if(GammaEnergy <= e2) { cs = pv->Value( 277 else if(GammaEnergy > e2) { cs = pv->Value( 278 279 return cs; 280 } 281 282 //******************************************** 283 284 void G4LowEPPolarizedComptonModel::SampleSecon 285 const G4MaterialCutsCouple* c 286 const G4DynamicParticle* aDyn 287 G4double, G4double) 288 { 289 290 //Determine number of digits (in decimal bas 291 G4double g4d_order = G4double(numeric_limits 292 G4double g4d_limit = std::pow(10.,-g4d_order 293 294 // The scattered gamma energy is sampled acc 295 // then accepted or rejected depending on th 296 // by factor from Klein - Nishina formula. 297 // Expression of the angular distribution as 298 // angular and energy distribution and Scatt 299 // D. E. Cullen "A simple model of photon tr 300 // Phys. Res. B 101 (1995). Method of sampli 301 // data are interpolated while in the articl 302 // Reference to the article is from J. Stepa 303 // and Electron Interaction Data for GEANT i 304 // TeV (draft). 305 // The random number techniques of Butcher & 306 // (Nucl Phys 20(1960),15). 307 308 G4double photonEnergy0 = aDynamicGamma->GetK 309 310 if (verboseLevel > 3) { 311 G4cout << "G4LowEPPolarizedComptonModel::S 312 << photonEnergy0/MeV << " in " << c 313 << G4endl; 314 } 315 // do nothing below the threshold 316 // should never get here because the XS is z 317 if (photonEnergy0 < LowEnergyLimit()) 318 return ; 319 320 G4double e0m = photonEnergy0 / electron_mass 321 G4ParticleMomentum photonDirection0 = aDynam 322 323 // Polarisation: check orientation of photon 324 // Fix if needed 325 326 G4ThreeVector photonPolarization0 = aDynamic 327 // Check if polarisation vector is perpendic 328 329 if (!(photonPolarization0.isOrthogonal(photo 330 { 331 photonPolarization0 = GetRandomPolarizat 332 } 333 else 334 { 335 if ((photonPolarization0.howOrthogonal(p 336 { 337 photonPolarization0 = GetPerpendicul 338 } 339 } 340 341 // Select randomly one element in the curren 342 const G4ParticleDefinition* particle = aDyn 343 const G4Element* elm = SelectRandomAtom(coup 344 G4int Z = (G4int)elm->GetZ(); 345 346 G4double LowEPPCepsilon0 = 1. / (1. + 2. * e 347 G4double LowEPPCepsilon0Sq = LowEPPCepsilon0 348 G4double alpha1 = -std::log(LowEPPCepsilon0) 349 G4double alpha2 = 0.5 * (1. - LowEPPCepsilon 350 351 G4double wlPhoton = h_Planck*c_light/photonE 352 353 // Sample the energy of the scattered photon 354 G4double LowEPPCepsilon; 355 G4double LowEPPCepsilonSq; 356 G4double oneCosT; 357 G4double sinT2; 358 G4double gReject; 359 360 if (verboseLevel > 3) { 361 G4cout << "Started loop to sample gamma en 362 } 363 364 do 365 { 366 if ( alpha1/(alpha1+alpha2) > G4UniformR 367 { 368 LowEPPCepsilon = G4Exp(-alpha1 * G4U 369 LowEPPCepsilonSq = LowEPPCepsilon * 370 } 371 else 372 { 373 LowEPPCepsilonSq = LowEPPCepsilon0Sq 374 LowEPPCepsilon = std::sqrt(LowEPPCep 375 } 376 377 oneCosT = (1. - LowEPPCepsilon) / ( LowE 378 sinT2 = oneCosT * (2. - oneCosT); 379 G4double x = std::sqrt(oneCosT/2.) / (wl 380 G4double scatteringFunction = ComputeSca 381 gReject = (1. - LowEPPCepsilon * sinT2 / 382 383 } while(gReject < G4UniformRand()*Z); 384 385 G4double cosTheta = 1. - oneCosT; 386 G4double sinTheta = std::sqrt(sinT2); 387 G4double phi = SetPhi(LowEPPCepsilon,sinT2); 388 G4double dirx = sinTheta * std::cos(phi); 389 G4double diry = sinTheta * std::sin(phi); 390 G4double dirz = cosTheta ; 391 392 // Set outgoing photon polarization 393 394 G4ThreeVector photonPolarization1 = SetNewPo 395 sinT2, 396 phi, 397 cosTheta); 398 399 // Scatter photon energy and Compton electro 400 // J. M. C. Brown, M. R. Dimmock, J. E. Gill 401 // "A low energy bound atomic electron Compt 402 // NIMB, Vol. 338, 77-88, 2014. 403 404 // Set constants and initialize scattering p 405 const G4double vel_c = c_light / (m/s); 406 const G4double momentum_au_to_nat = halfpi* 407 const G4double e_mass_kg = electron_mass_c2 408 409 const G4int maxDopplerIterations = 1000; 410 G4double bindingE = 0.; 411 G4double pEIncident = photonEnergy0 ; 412 G4double pERecoil = -1.; 413 G4double eERecoil = -1.; 414 G4double e_alpha =0.; 415 G4double e_beta = 0.; 416 417 G4double CE_emission_flag = 0.; 418 G4double ePAU = -1; 419 G4int shellIdx = 0; 420 G4double u_temp = 0; 421 G4double cosPhiE =0; 422 G4double sinThetaE =0; 423 G4double cosThetaE =0; 424 G4int iteration = 0; 425 426 if (verboseLevel > 3) { 427 G4cout << "Started loop to sample photon e 428 } 429 430 do{ 431 // *************************************** 432 // | Determine scatter photon energy 433 // *************************************** 434 do 435 { 436 iteration++; 437 438 // ***************************************** 439 // | Sample bound electron information 440 // ***************************************** 441 442 // Select shell based on shell occupancy 443 shellIdx = shellData->SelectRandomShell(Z); 444 bindingE = shellData->BindingEnergy(Z,shellI 445 446 // Randomly sample bound electron momentum ( 447 ePAU = profileData->RandomSelectMomentum(Z,s 448 449 // Convert to SI units 450 G4double ePSI = ePAU * momentum_au_to_nat; 451 452 //Calculate bound electron velocity and norm 453 u_temp = sqrt( ((ePSI*ePSI)*(vel_c*vel_c)) / 454 455 // Sample incident electron direction, amorp 456 e_alpha = pi*G4UniformRand(); 457 e_beta = twopi*G4UniformRand(); 458 459 // Total energy of system 460 G4double eEIncident = electron_mass_c2 / sqr 461 G4double systemE = eEIncident + pEIncident; 462 463 G4double gamma_temp = 1.0 / sqrt( 1 - (u_tem 464 G4double numerator = gamma_temp*electron_mas 465 G4double subdenom1 = u_temp*cosTheta*std::c 466 G4double subdenom2 = u_temp*sinTheta*std::si 467 G4double denominator = (1.0 - cosTheta) + ( 468 pERecoil = (numerator/denominator); 469 eERecoil = systemE - pERecoil; 470 CE_emission_flag = pEIncident - pERecoil; 471 } while ( (iteration <= maxDopplerIterat 472 473 // End of recalculation of photon energy w 474 // *************************************** 475 // | Determine ejected Compton electro 476 // *************************************** 477 478 // Calculate velocity of ejected Compton e 479 480 G4double a_temp = eERecoil / electron_mass 481 G4double u_p_temp = sqrt(1 - (1 / (a_temp* 482 483 // Coefficients and terms from simulatenou 484 G4double sinAlpha = std::sin(e_alpha); 485 G4double cosAlpha = std::cos(e_alpha); 486 G4double sinBeta = std::sin(e_beta); 487 G4double cosBeta = std::cos(e_beta); 488 489 G4double gamma = 1.0 / sqrt(1 - (u_temp*u_ 490 G4double gamma_p = 1.0 / sqrt(1 - (u_p_tem 491 492 G4double var_A = pERecoil*u_p_temp*sinThet 493 G4double var_B = u_p_temp* (pERecoil*cosTh 494 G4double var_C = (pERecoil-pEIncident) - ( 495 496 G4double var_D1 = gamma*electron_mass_c2*p 497 G4double var_D2 = (1 - (u_temp*cosTheta*co 498 G4double var_D3 = ((electron_mass_c2*elect 499 G4double var_D = var_D1*var_D2 + var_D3; 500 501 G4double var_E1 = ((gamma*gamma_p)*(electr 502 G4double var_E2 = gamma_p*electron_mass_c2 503 G4double var_E = var_E1 - var_E2; 504 505 G4double var_F1 = ((gamma*gamma_p)*(electr 506 G4double var_F2 = (gamma_p*electron_mass_c 507 G4double var_F = var_F1 - var_F2; 508 509 G4double var_G = (gamma*gamma_p)*(electron 510 511 // Two equations form a quadratic form of 512 // Coefficents and solution to quadratic 513 G4double var_W1 = (var_F*var_B - var_E*var 514 G4double var_W2 = (var_G*var_G)*(var_A*var 515 G4double var_W = var_W1 + var_W2; 516 517 G4double var_Y = 2.0*(((var_A*var_D-var_F* 518 519 G4double var_Z1 = (var_A*var_D - var_F*var 520 G4double var_Z2 = (var_G*var_G)*(var_C*var 521 G4double var_Z = var_Z1 + var_Z2; 522 G4double diff1 = var_Y*var_Y; 523 G4double diff2 = 4*var_W*var_Z; 524 G4double diff = diff1 - diff2; 525 526 // Check if diff is less than zero, if so 527 //Confirm that diff less than zero is due 528 //than 10^(-g4d_order), then set diff to z 529 530 if ((diff < 0.0) && (abs(diff / diff1) < g 531 { 532 diff = 0.0; 533 } 534 535 // Plus and minus of quadratic 536 G4double X_p = (-var_Y + sqrt (diff))/(2*v 537 G4double X_m = (-var_Y - sqrt (diff))/(2*v 538 539 // Floating point precision protection 540 // Check if X_p and X_m are greater than o 541 // Issue due to propagation of FPE and onl 542 if(X_p >1){X_p=1;} if(X_p<-1){X_p=-1;} 543 if(X_m >1){X_m=1;} if(X_m<-1){X_m=-1;} 544 // End of FP protection 545 546 G4double ThetaE = 0.; 547 548 // Randomly sample one of the two possible 549 G4double sol_select = G4UniformRand(); 550 551 if (sol_select < 0.5) 552 { 553 ThetaE = std::acos(X_p); 554 } 555 if (sol_select > 0.5) 556 { 557 ThetaE = std::acos(X_m); 558 } 559 560 cosThetaE = std::cos(ThetaE); 561 sinThetaE = std::sin(ThetaE); 562 G4double Theta = std::acos(cosTheta); 563 564 //Calculate electron Phi 565 G4double iSinThetaE = std::sqrt(1+std::tan 566 G4double iSinTheta = std::sqrt(1+std::tan( 567 G4double ivar_A = iSinTheta/ (pERecoil*u_p 568 // Trigs 569 cosPhiE = (var_C - var_B*cosThetaE)*(ivar_ 570 // End of calculation of ejection Compton 571 //Fix for floating point errors 572 } while ( (iteration <= maxDopplerIterations 573 574 // Revert to original if maximum number of i 575 if (iteration >= maxDopplerIterations) 576 { 577 pERecoil = photonEnergy0 ; 578 bindingE = 0.; 579 dirx=0.0; 580 diry=0.0; 581 dirz=1.0; 582 } 583 584 // Set "scattered" photon direction and ener 585 G4ThreeVector photonDirection1(dirx,diry,dir 586 SystemOfRefChange(photonDirection0,photonDir 587 photonPolarization0,photonPolarization 588 589 if (pERecoil > 0.) 590 { 591 fParticleChange->SetProposedKineticEnerg 592 fParticleChange->ProposeMomentumDirectio 593 fParticleChange->ProposePolarization(pho 594 595 // Set ejected Compton electron directio 596 G4double PhiE = std::acos(cosPhiE); 597 G4double eDirX = sinThetaE * std::cos(ph 598 G4double eDirY = sinThetaE * std::sin(ph 599 G4double eDirZ = cosThetaE; 600 601 G4double eKineticEnergy = pEIncident - p 602 603 G4ThreeVector eDirection(eDirX,eDirY,eDi 604 SystemOfRefChangeElect(photonDirection0, 605 photonPolarization0); 606 607 G4DynamicParticle* dp = new G4DynamicPar 608 eDirection,eKineticEnergy) ; 609 fvect->push_back(dp); 610 } 611 else 612 { 613 fParticleChange->SetProposedKineticEnerg 614 fParticleChange->ProposeTrackStatus(fSto 615 } 616 617 // sample deexcitation 618 // 619 if (verboseLevel > 3) { 620 G4cout << "Started atomic de-excitation " 621 } 622 623 if(fAtomDeexcitation && iteration < maxDoppl 624 G4int index = couple->GetIndex(); 625 if(fAtomDeexcitation->CheckDeexcitationAct 626 std::size_t nbefore = fvect->size(); 627 G4AtomicShellEnumerator as = G4AtomicShe 628 const G4AtomicShell* shell = fAtomDeexci 629 fAtomDeexcitation->GenerateParticles(fve 630 std::size_t nafter = fvect->size(); 631 if(nafter > nbefore) { 632 for (std::size_t i=nbefore; i<nafter; 633 //Check if there is enough residual 634 if (bindingE >= ((*fvect)[i])->GetKi 635 { 636 //Ok, this is a valid secondary: keep 637 bindingE -= ((*fvect)[i])->GetKineticE 638 } 639 else 640 { 641 //Invalid secondary: not enough energy 642 //Keep its energy in the local deposit 643 delete (*fvect)[i]; 644 (*fvect)[i]=nullptr; 645 } 646 } 647 } 648 } 649 } 650 651 //This should never happen 652 if(bindingE < 0.0) 653 G4Exception("G4LowEPPolarizedComptonModel: 654 "em2051",FatalException,"Negative local en 655 656 fParticleChange->ProposeLocalEnergyDeposit(b 657 } 658 659 //******************************************** 660 661 G4double 662 G4LowEPPolarizedComptonModel::ComputeScatterin 663 { 664 G4double value = Z; 665 if (x <= ScatFuncFitParam[Z][2]) { 666 667 G4double lgq = G4Log(x)/ln10; 668 669 if (lgq < ScatFuncFitParam[Z][1]) { 670 value = ScatFuncFitParam[Z][3] + lgq*Sca 671 } else { 672 value = ScatFuncFitParam[Z][5] + lgq*Sca 673 lgq*lgq*ScatFuncFitParam[Z][7] + lgq*l 674 } 675 value = G4Exp(value*ln10); 676 } 677 return value; 678 } 679 680 681 //******************************************** 682 683 void 684 G4LowEPPolarizedComptonModel::InitialiseForEle 685 G4int Z) 686 { 687 G4AutoLock l(&LowEPPolarizedComptonModelMute 688 if(!data[Z]) { ReadData(Z); } 689 l.unlock(); 690 } 691 692 //******************************************** 693 694 //Fitting data to compute scattering function 695 const G4double G4LowEPPolarizedComptonModel::S 696 { 0, 0., 0., 0., 0., 697 { 1, 6.673, 1.49968E+08, -14.352, 1.999, -1 698 { 2, 6.500, 2.50035E+08, -14.215, 1.970, -5 699 { 3, 6.551, 3.99945E+08, -13.555, 1.993, -6 700 { 4, 6.500, 5.00035E+08, -13.746, 1.998, -1 701 { 5, 6.500, 5.99791E+08, -13.800, 1.998, -1 702 { 6, 6.708, 6.99842E+08, -13.885, 1.999, -1 703 { 7, 6.685, 7.99834E+08, -13.885, 2.000, -1 704 { 8, 6.669, 7.99834E+08, -13.962, 2.001, -1 705 { 9, 6.711, 7.99834E+08, -13.999, 2.000, -1 706 { 10, 6.702, 7.99834E+08, -14.044, 1.999, -1 707 { 11, 6.425, 1.00000E+09, -13.423, 1.993, -4 708 { 12, 6.542, 1.00000E+09, -13.389, 1.997, -5 709 { 13, 6.570, 1.49968E+09, -13.401, 1.997, -6 710 { 14, 6.364, 1.49968E+09, -13.452, 1.999, -7 711 { 15, 6.500, 1.49968E+09, -13.488, 1.998, -8 712 { 16, 6.500, 1.49968E+09, -13.532, 1.998, -9 713 { 17, 6.500, 1.49968E+09, -13.584, 2.000, -9 714 { 18, 6.500, 1.49968E+09, -13.618, 1.999, -1 715 { 19, 6.500, 1.99986E+09, -13.185, 1.992, -5 716 { 20, 6.490, 1.99986E+09, -13.123, 1.993, -5 717 { 21, 6.498, 1.99986E+09, -13.157, 1.994, -5 718 { 22, 6.495, 1.99986E+09, -13.183, 1.994, -5 719 { 23, 6.487, 1.99986E+09, -13.216, 1.995, -5 720 { 24, 6.500, 1.99986E+09, -13.330, 1.997, -6 721 { 25, 6.488, 1.99986E+09, -13.277, 1.997, -5 722 { 26, 6.500, 5.00035E+09, -13.292, 1.997, -6 723 { 27, 6.500, 5.00035E+09, -13.321, 1.998, -6 724 { 28, 6.500, 5.00035E+09, -13.340, 1.998, -6 725 { 29, 6.500, 5.00035E+09, -13.439, 1.998, -6 726 { 30, 6.500, 5.00035E+09, -13.383, 1.999, -6 727 { 31, 6.500, 5.00035E+09, -13.349, 1.997, -6 728 { 32, 6.500, 5.00035E+09, -13.373, 1.999, -6 729 { 33, 6.500, 5.00035E+09, -13.395, 1.999, -6 730 { 34, 6.500, 5.00035E+09, -13.417, 1.999, -6 731 { 35, 6.500, 5.00035E+09, -13.442, 2.000, -7 732 { 36, 6.500, 5.00035E+09, -13.451, 1.998, -7 733 { 37, 6.500, 5.00035E+09, -13.082, 1.991, -4 734 { 38, 6.465, 5.00035E+09, -13.022, 1.993, -4 735 { 39, 6.492, 5.00035E+09, -13.043, 1.994, -4 736 { 40, 6.499, 5.00035E+09, -13.064, 1.994, -4 737 { 41, 6.384, 5.00035E+09, -13.156, 1.996, -5 738 { 42, 6.500, 5.00035E+09, -13.176, 1.996, -5 739 { 43, 6.500, 5.00035E+09, -13.133, 1.997, -5 740 { 44, 6.500, 5.00035E+09, -13.220, 1.996, -5 741 { 45, 6.500, 5.00035E+09, -13.246, 1.998, -5 742 { 46, 6.500, 5.00035E+09, -13.407, 1.999, -7 743 { 47, 6.500, 5.00035E+09, -13.277, 1.998, -6 744 { 48, 6.500, 5.00035E+09, -13.222, 1.998, -5 745 { 49, 6.500, 5.00035E+09, -13.199, 1.997, -5 746 { 50, 6.500, 5.00035E+09, -13.215, 1.998, -5 747 { 51, 6.500, 5.00035E+09, -13.230, 1.998, -6 748 { 52, 6.500, 7.99834E+09, -13.246, 1.998, -6 749 { 53, 6.500, 5.00035E+09, -13.262, 1.998, -6 750 { 54, 6.500, 7.99834E+09, -13.279, 1.998, -6 751 { 55, 6.500, 5.00035E+09, -12.951, 1.990, -4 752 { 56, 6.425, 5.00035E+09, -12.882, 1.992, -3 753 { 57, 6.466, 2.82488E+09, -12.903, 1.992, -3 754 { 58, 6.451, 5.00035E+09, -12.915, 1.993, -4 755 { 59, 6.434, 5.00035E+09, -12.914, 1.993, -4 756 { 60, 6.444, 5.00035E+09, -12.922, 1.992, -3 757 { 61, 6.414, 7.99834E+09, -12.930, 1.993, -4 758 { 62, 6.420, 7.99834E+09, -12.938, 1.992, -4 759 { 63, 6.416, 7.99834E+09, -12.946, 1.993, -4 760 { 64, 6.443, 7.99834E+09, -12.963, 1.993, -4 761 { 65, 6.449, 7.99834E+09, -12.973, 1.993, -4 762 { 66, 6.419, 7.99834E+09, -12.966, 1.993, -4 763 { 67, 6.406, 7.99834E+09, -12.976, 1.993, -4 764 { 68, 6.424, 7.99834E+09, -12.986, 1.993, -4 765 { 69, 6.417, 7.99834E+09, -12.989, 1.993, -4 766 { 70, 6.405, 7.99834E+09, -13.000, 1.994, -4 767 { 71, 6.449, 7.99834E+09, -13.015, 1.994, -4 768 { 72, 6.465, 7.99834E+09, -13.030, 1.994, -4 769 { 73, 6.447, 7.99834E+09, -13.048, 1.996, -4 770 { 74, 6.452, 7.99834E+09, -13.073, 1.997, -4 771 { 75, 6.432, 7.99834E+09, -13.082, 1.997, -4 772 { 76, 6.439, 7.99834E+09, -13.100, 1.997, -4 773 { 77, 6.432, 7.99834E+09, -13.110, 1.997, -4 774 { 78, 6.500, 7.99834E+09, -13.185, 1.997, -5 775 { 79, 6.500, 7.99834E+09, -13.200, 1.997, -5 776 { 80, 6.500, 7.99834E+09, -13.156, 1.998, -4 777 { 81, 6.500, 7.99834E+09, -13.128, 1.997, -4 778 { 82, 6.500, 7.99834E+09, -13.134, 1.997, -5 779 { 83, 6.500, 7.99834E+09, -13.148, 1.998, -5 780 { 84, 6.500, 7.99834E+09, -13.161, 1.998, -5 781 { 85, 6.500, 7.99834E+09, -13.175, 1.998, -5 782 { 86, 6.500, 7.99834E+09, -13.189, 1.998, -5 783 { 87, 6.500, 7.99834E+09, -12.885, 1.990, -3 784 { 88, 6.417, 7.99834E+09, -12.816, 1.991, -3 785 { 89, 6.442, 7.99834E+09, -12.831, 1.992, -3 786 { 90, 6.463, 7.99834E+09, -12.850, 1.993, -3 787 { 91, 6.447, 7.99834E+09, -12.852, 1.993, -3 788 { 92, 6.439, 7.99834E+09, -12.858, 1.993, -3 789 { 93, 6.437, 1.00000E+10, -12.866, 1.993, -3 790 { 94, 6.432, 7.99834E+09, -12.862, 1.993, -3 791 { 95, 6.435, 7.99834E+09, -12.869, 1.993, -3 792 { 96, 6.449, 1.00000E+10, -12.886, 1.993, -3 793 { 97, 6.446, 1.00000E+10, -12.892, 1.993, -4 794 { 98, 6.421, 1.00000E+10, -12.887, 1.993, -3 795 { 99, 6.414, 1.00000E+10, -12.894, 1.993, -3 796 {100, 6.412, 1.00000E+10, -12.900, 1.993, -3 797 }; 798 799 //******************************************** 800 801 //Supporting functions for photon polarisation 802 G4double G4LowEPPolarizedComptonModel::SetPhi( 803 G4double sinT2) 804 { 805 G4double rand1; 806 G4double rand2; 807 G4double phiProbability; 808 G4double phi; 809 G4double a, b; 810 811 do 812 { 813 rand1 = G4UniformRand(); 814 rand2 = G4UniformRand(); 815 phiProbability=0.; 816 phi = twopi*rand1; 817 818 a = 2*sinT2; 819 b = energyRate + 1/energyRate; 820 821 phiProbability = 1 - (a/b)*(std::cos(phi 822 } 823 while ( rand2 > phiProbability ); 824 return phi; 825 } 826 827 //******************************************** 828 829 G4ThreeVector G4LowEPPolarizedComptonModel::Se 830 { 831 G4double dx = a.x(); 832 G4double dy = a.y(); 833 G4double dz = a.z(); 834 G4double x = dx < 0.0 ? -dx : dx; 835 G4double y = dy < 0.0 ? -dy : dy; 836 G4double z = dz < 0.0 ? -dz : dz; 837 if (x < y) { 838 return x < z ? G4ThreeVector(-dy,dx,0) : G 839 }else{ 840 return y < z ? G4ThreeVector(dz,0,-dx) : G 841 } 842 } 843 844 //******************************************** 845 846 G4ThreeVector G4LowEPPolarizedComptonModel::Ge 847 { 848 G4ThreeVector d0 = direction0.unit(); 849 G4ThreeVector a1 = SetPerpendicularVector(d0 850 G4ThreeVector a0 = a1.unit(); // unit vector 851 852 G4double rand1 = G4UniformRand(); 853 854 G4double angle = twopi*rand1; // random pola 855 G4ThreeVector b0 = d0.cross(a0); // cross pr 856 857 G4ThreeVector c; 858 859 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 860 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 861 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 862 863 G4ThreeVector c0 = c.unit(); 864 865 return c0; 866 } 867 868 //******************************************** 869 870 G4ThreeVector G4LowEPPolarizedComptonModel::Ge 871 (const G4ThreeVector& photonDirection, const G 872 { 873 // 874 // The polarization of a photon is always pe 875 // Therefore this function removes those vec 876 // points in direction of photonDirection 877 // 878 // Mathematically we search the projection o 879 // plains normal vector. 880 // The basic equation can be found in each g 881 // p = a - (a o n)/(n o n)*n 882 883 return photonPolarization - photonPolarizati 884 } 885 886 //******************************************** 887 888 G4ThreeVector G4LowEPPolarizedComptonModel::Se 889 G4double sinT2, 890 G4double phi, 891 G4double costheta) 892 { 893 G4double rand1; 894 G4double rand2; 895 G4double cosPhi = std::cos(phi); 896 G4double sinPhi = std::sin(phi); 897 G4double sinTheta = std::sqrt(sinT2); 898 G4double cosP2 = cosPhi*cosPhi; 899 G4double normalisation = std::sqrt(1. - cosP 900 901 // Method based on: 902 // D. Xu, Z. He and F. Zhang 903 // "Detection of Gamma Ray Polarization Usin 904 // IEEE TNS, Vol. 52(4), 1160-1164, 2005. 905 906 // Determination of Theta 907 908 G4double theta; 909 rand1 = G4UniformRand(); 910 rand2 = G4UniformRand(); 911 912 if (rand1<(LowEPPCepsilon+1.0/LowEPPCepsilon 913 { 914 if (rand2<0.5) 915 theta = pi/2.0; 916 else 917 theta = 3.0*pi/2.0; 918 } 919 else 920 { 921 if (rand2<0.5) 922 theta = 0; 923 else 924 theta = pi; 925 } 926 G4double cosBeta = std::cos(theta); 927 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 928 929 G4ThreeVector photonPolarization1; 930 931 G4double xParallel = normalisation*cosBeta; 932 G4double yParallel = -(sinT2*cosPhi*sinPhi)* 933 G4double zParallel = -(costheta*sinTheta*cos 934 G4double xPerpendicular = 0.; 935 G4double yPerpendicular = (costheta)*sinBeta 936 G4double zPerpendicular = -(sinTheta*sinPhi) 937 938 G4double xTotal = (xParallel + xPerpendicula 939 G4double yTotal = (yParallel + yPerpendicula 940 G4double zTotal = (zParallel + zPerpendicula 941 942 photonPolarization1.setX(xTotal); 943 photonPolarization1.setY(yTotal); 944 photonPolarization1.setZ(zTotal); 945 946 return photonPolarization1; 947 948 } 949 950 //******************************************** 951 void G4LowEPPolarizedComptonModel::SystemOfRef 952 G4ThreeVector& direction1, 953 G4ThreeVector& polarization0, 954 G4ThreeVector& polarization1) 955 { 956 // direction0 is the original photon directi 957 // polarization0 is the original photon pola 958 // need to specify y axis in the real refere 959 G4ThreeVector Axis_Z0 = direction0.unit(); 960 G4ThreeVector Axis_X0 = polarization0.unit() 961 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 962 963 G4double direction_x = direction1.getX(); 964 G4double direction_y = direction1.getY(); 965 G4double direction_z = direction1.getZ(); 966 967 direction1 = (direction_x*Axis_X0 + directio 968 G4double polarization_x = polarization1.getX 969 G4double polarization_y = polarization1.getY 970 G4double polarization_z = polarization1.getZ 971 972 polarization1 = (polarization_x*Axis_X0 + po 973 974 } 975 976 //******************************************** 977 void G4LowEPPolarizedComptonModel::SystemOfRef 978 G4ThreeVector& edirection, 979 G4ThreeVector& ppolarization) 980 { 981 // direction0 is the original photon directi 982 // polarization0 is the original photon pola 983 // need to specify y axis in the real refere 984 G4ThreeVector Axis_Z0 = pdirection.unit(); 985 G4ThreeVector Axis_X0 = ppolarization.unit() 986 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 987 988 G4double direction_x = edirection.getX(); 989 G4double direction_y = edirection.getY(); 990 G4double direction_z = edirection.getZ(); 991 992 edirection = (direction_x*Axis_X0 + directio 993 } 994 995 996