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 // Author: Luciano Pandola 28 // 29 // History: 30 // -------- 31 // 08 Jan 2010 L Pandola First implementati 32 // 01 Feb 2011 L Pandola Suppress fake ener 33 // is active. 34 // Make sure that flu 35 // only if above thre 36 // 25 May 2011 L Pandola Renamed (make v200 37 // 10 Jun 2011 L Pandola Migrate atomic dee 38 // 07 Oct 2011 L Pandola Bug fix (potential 39 // 27 Sep 2013 L Pandola Migrate to MT para 40 // tables. 41 // 02 Oct 2013 L Pandola Rewrite sampling a 42 // to improve CPU per 43 // 44 45 #include "G4PenelopePhotoElectricModel.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4ParticleDefinition.hh" 49 #include "G4MaterialCutsCouple.hh" 50 #include "G4DynamicParticle.hh" 51 #include "G4PhysicsTable.hh" 52 #include "G4PhysicsFreeVector.hh" 53 #include "G4ElementTable.hh" 54 #include "G4Element.hh" 55 #include "G4AtomicTransitionManager.hh" 56 #include "G4AtomicShell.hh" 57 #include "G4Gamma.hh" 58 #include "G4Electron.hh" 59 #include "G4AutoLock.hh" 60 #include "G4LossTableManager.hh" 61 #include "G4Exp.hh" 62 63 //....oooOO0OOooo........oooOO0OOooo........oo 64 65 const G4int G4PenelopePhotoElectricModel::fMax 66 G4PhysicsTable* G4PenelopePhotoElectricModel:: 67 68 //....oooOO0OOooo........oooOO0OOooo........oo 69 70 G4PenelopePhotoElectricModel::G4PenelopePhotoE 71 const G4String& nam) 72 :G4VEmModel(nam),fParticleChange(nullptr),fP 73 fAtomDeexcitation(nullptr),fIsInitialised(f 74 { 75 fIntrinsicLowEnergyLimit = 100.0*eV; 76 fIntrinsicHighEnergyLimit = 100.0*GeV; 77 // SetLowEnergyLimit(fIntrinsicLowEnergyLim 78 SetHighEnergyLimit(fIntrinsicHighEnergyLimit 79 // 80 81 if (part) 82 SetParticle(part); 83 84 fVerboseLevel= 0; 85 // Verbosity scale: 86 // 0 = nothing 87 // 1 = warning for energy non-conservation 88 // 2 = details of energy budget 89 // 3 = calculation of cross sections, file o 90 // 4 = entering in methods 91 92 //Mark this model as "applicable" for atomic 93 SetDeexcitationFlag(true); 94 95 fTransitionManager = G4AtomicTransitionManag 96 } 97 98 //....oooOO0OOooo........oooOO0OOooo........oo 99 100 G4PenelopePhotoElectricModel::~G4PenelopePhoto 101 { 102 if (IsMaster() || fLocalTable) 103 { 104 for(G4int i=0; i<=fMaxZ; ++i) 105 { 106 if(fLogAtomicShellXS[i]) { 107 fLogAtomicShellXS[i]->clearAndDestroy(); 108 delete fLogAtomicShellXS[i]; 109 fLogAtomicShellXS[i] = nullptr; 110 } 111 } 112 } 113 } 114 115 //....oooOO0OOooo........oooOO0OOooo........oo 116 117 void G4PenelopePhotoElectricModel::Initialise( 118 const G4DataVector& cuts) 119 { 120 if (fVerboseLevel > 3) 121 G4cout << "Calling G4PenelopePhotoElectri 122 123 fAtomDeexcitation = G4LossTableManager::Inst 124 //Issue warning if the AtomicDeexcitation ha 125 if (!fAtomDeexcitation) 126 { 127 G4cout << G4endl; 128 G4cout << "WARNING from G4PenelopePhotoE 129 G4cout << "Atomic de-excitation module i 130 G4cout << "any fluorescence/Auger emissi 131 G4cout << "Please make sure this is inte 132 } 133 134 SetParticle(particle); 135 136 //Only the master model creates/fills/destro 137 if (IsMaster() && particle == fParticle) 138 { 139 G4ProductionCutsTable* theCoupleTable = 140 G4ProductionCutsTable::GetProductionCutsTabl 141 142 for (G4int i=0;i<(G4int)theCoupleTable-> 143 { 144 const G4Material* material = 145 theCoupleTable->GetMaterialCutsCouple(i) 146 const G4ElementVector* theElementVector = 147 148 for (std::size_t j=0;j<material->GetNumber 149 { 150 G4int iZ = theElementVector->at(j)->Ge 151 //read data files only in the master 152 if (!fLogAtomicShellXS[iZ]) 153 ReadDataFile(iZ); 154 } 155 } 156 157 InitialiseElementSelectors(particle,cuts 158 159 if (fVerboseLevel > 0) { 160 G4cout << "Penelope Photo-Electric model v20 161 << "Energy range: " 162 << LowEnergyLimit() / MeV << " MeV - 163 << HighEnergyLimit() / GeV << " GeV"; 164 } 165 } 166 167 if(fIsInitialised) return; 168 fParticleChange = GetParticleChangeForGamma( 169 fIsInitialised = true; 170 171 } 172 173 void G4PenelopePhotoElectricModel::InitialiseL 174 G4VEmModel *masterModel) 175 { 176 if (fVerboseLevel > 3) 177 G4cout << "Calling G4PenelopePhotoElectri 178 // 179 //Check that particle matches: one might hav 180 //for e+ and e-). 181 // 182 if (part == fParticle) 183 { 184 SetElementSelectors(masterModel->GetElem 185 186 //Get the const table pointers from the 187 const G4PenelopePhotoElectricModel* theM 188 static_cast<G4PenelopePhotoElectricModel*> ( 189 for(G4int i=0; i<=fMaxZ; ++i) 190 fLogAtomicShellXS[i] = theModel->fLogAtomicS 191 //Same verbosity for all workers, as the 192 fVerboseLevel = theModel->fVerboseLevel; 193 } 194 195 return; 196 } 197 198 //....oooOO0OOooo........oooOO0OOooo........oo 199 namespace { G4Mutex PenelopePhotoElectricMode 200 G4double G4PenelopePhotoElectricModel::Compute 201 const G4ParticleDefinition*, 202 G4double energy, 203 G4double Z, G4double, 204 G4double, G4double) 205 { 206 // 207 // Penelope model v2008 208 // 209 if (fVerboseLevel > 3) 210 G4cout << "Calling ComputeCrossSectionPerA 211 212 G4int iZ = G4int(Z); 213 214 if (!fLogAtomicShellXS[iZ]) 215 { 216 //If we are here, it means that Initiali 217 //not filled up. This can happen in a Un 218 if (fVerboseLevel > 0) 219 { 220 //Issue a G4Exception (warning) only in ve 221 G4ExceptionDescription ed; 222 ed << "Unable to retrieve the shell cross 223 ed << "This can happen only in Unit Tests 224 G4Exception("G4PenelopePhotoElectricModel: 225 "em2038",JustWarning,ed); 226 } 227 //protect file reading via autolock 228 G4AutoLock lock(&PenelopePhotoElectricMo 229 ReadDataFile(iZ); 230 lock.unlock(); 231 } 232 233 G4double cross = 0; 234 G4PhysicsTable* theTable = fLogAtomicShellX 235 G4PhysicsFreeVector* totalXSLog = (G4Physics 236 237 if (!totalXSLog) 238 { 239 G4Exception("G4PenelopePhotoElectricMod 240 "em2039",FatalException, 241 "Unable to retrieve the total cross sec 242 return 0; 243 } 244 G4double logene = G4Log(energy); 245 G4double logXS = totalXSLog->Value(logene); 246 cross = G4Exp(logXS); 247 248 if (fVerboseLevel > 2) 249 G4cout << "Photoelectric cross section at 250 " = " << cross/barn << " barn" << G4endl 251 return cross; 252 } 253 254 //....oooOO0OOooo........oooOO0OOooo........oo 255 256 void G4PenelopePhotoElectricModel::SampleSecon 257 const G4MaterialCutsCouple* c 258 const G4DynamicParticle* aDyn 259 G4double, 260 G4double) 261 { 262 // 263 // Photoelectric effect, Penelope model v200 264 // 265 // The target atom and the target shell are 266 // database 267 // D.E. Cullen et al., Report UCRL-50400 (1 268 // The angular distribution of the electron 269 // according to the Sauter distribution from 270 // F. Sauter, Ann. Phys. 11 (1931) 454 271 // The energy of the final electron is given 272 // the binding energy. Fluorescence de-excit 273 // (to fill the vacancy) according to the ge 274 // J. Stepanek, Comp. Phys. Comm. 1206 pp 1 275 276 if (fVerboseLevel > 3) 277 G4cout << "Calling SamplingSecondaries() o 278 279 G4double photonEnergy = aDynamicGamma->GetKi 280 281 // always kill primary 282 fParticleChange->ProposeTrackStatus(fStopAnd 283 fParticleChange->SetProposedKineticEnergy(0. 284 285 if (photonEnergy <= fIntrinsicLowEnergyLimit 286 { 287 fParticleChange->ProposeLocalEnergyDepos 288 return ; 289 } 290 291 G4ParticleMomentum photonDirection = aDynami 292 293 // Select randomly one element in the curren 294 if (fVerboseLevel > 2) 295 G4cout << "Going to select element in " << 296 297 // atom can be selected efficiently if eleme 298 const G4Element* anElement = 299 SelectRandomAtom(couple,G4Gamma::GammaDefi 300 G4int Z = anElement->GetZasInt(); 301 if (fVerboseLevel > 2) 302 G4cout << "Selected " << anElement->GetNam 303 304 // Select the ionised shell in the current a 305 //shellIndex = 0 --> K shell 306 // 1-3 --> L shells 307 // 4-8 --> M shells 308 // 9 --> outer shells cumulative 309 // 310 std::size_t shellIndex = SelectRandomShell(Z 311 312 if (fVerboseLevel > 2) 313 G4cout << "Selected shell " << shellIndex 314 315 // Retrieve the corresponding identifier and 316 const G4AtomicTransitionManager* transitionM 317 318 //The number of shell cross section possibly 319 //might be different from the number of shel 320 //(namely, Penelope may contain more shell, 321 //In order to avoid a warning message from t 322 //add this protection. Results are anyway ch 323 //has a shellID>maxID, it sets the shellID t 324 std::size_t numberOfShells = (std::size_t) t 325 if (shellIndex >= numberOfShells) 326 shellIndex = numberOfShells-1; 327 328 const G4AtomicShell* shell = fTransitionMana 329 G4double bindingEnergy = shell->BindingEnerg 330 331 //Penelope considers only K, L and M shells. 332 //not included in the Penelope database. If 333 //shellIndex = 9, it means that an outer she 334 //Penelope recipe is to set bindingEnergy = 335 //to the electron) and to disregard fluoresc 336 if (shellIndex == 9) 337 bindingEnergy = 0.*eV; 338 339 G4double localEnergyDeposit = 0.0; 340 G4double cosTheta = 1.0; 341 342 // Primary outcoming electron 343 G4double eKineticEnergy = photonEnergy - bin 344 345 // There may be cases where the binding ener 346 // In such cases do not generate secondaries 347 if (eKineticEnergy > 0.) 348 { 349 // The electron is created 350 // Direction sampled from the Sauter dis 351 cosTheta = SampleElectronDirection(eKine 352 G4double sinTheta = std::sqrt(1-cosTheta 353 G4double phi = twopi * G4UniformRand() ; 354 G4double dirx = sinTheta * std::cos(phi) 355 G4double diry = sinTheta * std::sin(phi) 356 G4double dirz = cosTheta ; 357 G4ThreeVector electronDirection(dirx,dir 358 electronDirection.rotateUz(photonDirecti 359 G4DynamicParticle* electron = new G4Dyna 360 electronDirection, 361 eKineticEnergy); 362 fvect->push_back(electron); 363 } 364 else 365 bindingEnergy = photonEnergy; 366 367 G4double energyInFluorescence = 0; //testing 368 G4double energyInAuger = 0; //testing purpos 369 370 //Now, take care of fluorescence, if require 371 //recipe, I have to skip fluoresence complet 372 //(= sampling of a shell outer than K,L,M) 373 if (fAtomDeexcitation && shellIndex<9) 374 { 375 G4int index = couple->GetIndex(); 376 if (fAtomDeexcitation->CheckDeexcitation 377 { 378 std::size_t nBefore = fvect->size(); 379 fAtomDeexcitation->GenerateParticles(fvect 380 std::size_t nAfter = fvect->size(); 381 382 if (nAfter > nBefore) //actual production 383 { 384 for (std::size_t j=nBefore;j<nAfter;++ 385 { 386 G4double itsEnergy = ((*fvect)[j])->GetK 387 if (itsEnergy < bindingEnergy) // valid 388 { 389 bindingEnergy -= itsEnergy; 390 if (((*fvect)[j])->GetParticleDefini 391 energyInFluorescence += itsEnergy; 392 else if (((*fvect)[j])->GetParticleD 393 energyInAuger += itsEnergy; 394 } 395 else //invalid secondary: takes more tha 396 { 397 delete (*fvect)[j]; 398 (*fvect)[j] = nullptr; 399 } 400 } 401 } 402 } 403 } 404 405 //Residual energy is deposited locally 406 localEnergyDeposit += bindingEnergy; 407 408 if (localEnergyDeposit < 0) //Should not be: 409 { 410 G4Exception("G4PenelopePhotoElectricMode 411 "em2099",JustWarning,"WARNING: Negative 412 localEnergyDeposit = 0; 413 } 414 415 fParticleChange->ProposeLocalEnergyDeposit(l 416 417 if (fVerboseLevel > 1) 418 { 419 G4cout << "----------------------------- 420 G4cout << "Energy balance from G4Penelop 421 G4cout << "Selected shell: " << WriteTar 422 anElement->GetName() << G4endl; 423 G4cout << "Incoming photon energy: " << 424 G4cout << "----------------------------- 425 if (eKineticEnergy) 426 G4cout << "Outgoing electron " << eKineticEn 427 if (energyInFluorescence) 428 G4cout << "Fluorescence x-rays: " << energyI 429 if (energyInAuger) 430 G4cout << "Auger electrons: " << energyInAug 431 G4cout << "Local energy deposit " << loc 432 G4cout << "Total final state: " << 433 (eKineticEnergy+energyInFluorescence+localEn 434 " keV" << G4endl; 435 G4cout << "----------------------------- 436 } 437 if (fVerboseLevel > 0) 438 { 439 G4double energyDiff = 440 std::fabs(eKineticEnergy+energyInFluorescenc 441 if (energyDiff > 0.05*keV) 442 { 443 G4cout << "Warning from G4PenelopePhotoEle 444 (eKineticEnergy+energyInFluorescence+loc 445 << " keV (final) vs. " << 446 photonEnergy/keV << " keV (initial)" << 447 G4cout << "------------------------------- 448 G4cout << "Energy balance from G4PenelopeP 449 G4cout << "Selected shell: " << WriteTarge 450 anElement->GetName() << G4endl; 451 G4cout << "Incoming photon energy: " << ph 452 G4cout << "------------------------------- 453 if (eKineticEnergy) 454 G4cout << "Outgoing electron " << eKinet 455 if (energyInFluorescence) 456 G4cout << "Fluorescence x-rays: " << ene 457 if (energyInAuger) 458 G4cout << "Auger electrons: " << energyI 459 G4cout << "Local energy deposit " << local 460 G4cout << "Total final state: " << 461 (eKineticEnergy+energyInFluorescence+loc 462 " keV" << G4endl; 463 G4cout << "------------------------------- 464 } 465 } 466 } 467 468 //....oooOO0OOooo........oooOO0OOooo........oo 469 470 G4double G4PenelopePhotoElectricModel::SampleE 471 { 472 G4double costheta = 1.0; 473 if (energy>1*GeV) return costheta; 474 475 //1) initialize energy-dependent variables 476 // Variable naming according to Eq. (2.24) o 477 // (pag. 44) 478 G4double gamma = 1.0 + energy/electron_mass_ 479 G4double gamma2 = gamma*gamma; 480 G4double beta = std::sqrt((gamma2-1.0)/gamma 481 482 // ac corresponds to "A" of Eq. (2.31) 483 // 484 G4double ac = (1.0/beta) - 1.0; 485 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(ga 486 G4double a2 = ac + 2.0; 487 G4double gtmax = 2.0*(a1 + 1.0/ac); 488 489 G4double tsam = 0; 490 G4double gtr = 0; 491 492 //2) sampling. Eq. (2.31) of Penelope Manual 493 // tsam = 1-std::cos(theta) 494 // gtr = rejection function according to Eq. 495 do{ 496 G4double rand = G4UniformRand(); 497 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(r 498 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam)); 499 }while(G4UniformRand()*gtmax > gtr); 500 costheta = 1.0-tsam; 501 502 return costheta; 503 } 504 505 //....oooOO0OOooo........oooOO0OOooo........oo 506 507 void G4PenelopePhotoElectricModel::ReadDataFil 508 { 509 if (!IsMaster()) 510 //Should not be here! 511 G4Exception("G4PenelopePhotoElectricModel: 512 "em0100",FatalException,"Worker thread in 513 514 if (fVerboseLevel > 2) 515 { 516 G4cout << "G4PenelopePhotoElectricModel: 517 G4cout << "Going to read PhotoElectric d 518 } 519 520 const char* path = G4FindDataDir("G4LEDATA 521 if(!path) 522 { 523 G4String excep = "G4PenelopePhotoElectri 524 G4Exception("G4PenelopePhotoElectricMode 525 "em0006",FatalException,excep); 526 return; 527 } 528 529 /* 530 Read the cross section file 531 */ 532 std::ostringstream ost; 533 if (Z>9) 534 ost << path << "/penelope/photoelectric/pd 535 else 536 ost << path << "/penelope/photoelectric/pd 537 std::ifstream file(ost.str().c_str()); 538 if (!file.is_open()) 539 { 540 G4String excep = "G4PenelopePhotoElectri 541 G4Exception("G4PenelopePhotoElectricMode 542 "em0003",FatalException,excep); 543 } 544 //I have to know in advance how many points 545 //to initialize the G4PhysicsFreeVector() 546 std::size_t ndata=0; 547 G4String line; 548 while( getline(file, line) ) 549 ndata++; 550 ndata -= 1; 551 //G4cout << "Found: " << ndata << " lines" < 552 553 file.clear(); 554 file.close(); 555 file.open(ost.str().c_str()); 556 557 G4int readZ =0; 558 std::size_t nShells= 0; 559 file >> readZ >> nShells; 560 561 if (fVerboseLevel > 3) 562 G4cout << "Element Z=" << Z << " , nShells 563 564 //check the right file is opened. 565 if (readZ != Z || nShells <= 0 || nShells > 566 { 567 G4ExceptionDescription ed; 568 ed << "Corrupted data file for Z=" << Z 569 G4Exception("G4PenelopePhotoElectricMode 570 "em0005",FatalException,ed); 571 return; 572 } 573 G4PhysicsTable* thePhysicsTable = new G4Phys 574 575 //the table has to contain nShell+1 G4Physic 576 //(theTable)[0] --> total cross section 577 //(theTable)[ishell] --> cross section for s 578 579 //reserve space for the vectors 580 //everything is log-log 581 for (std::size_t i=0;i<nShells+1;++i) 582 thePhysicsTable->push_back(new G4PhysicsFr 583 584 std::size_t k =0; 585 for (k=0;k<ndata && !file.eof();++k) 586 { 587 G4double energy = 0; 588 G4double aValue = 0; 589 file >> energy ; 590 energy *= eV; 591 G4double logene = G4Log(energy); 592 //loop on the columns 593 for (std::size_t i=0;i<nShells+1;++i) 594 { 595 file >> aValue; 596 aValue *= barn; 597 G4PhysicsFreeVector* theVec = (G4PhysicsFr 598 if (aValue < 1e-40*cm2) //protection again 599 aValue = 1e-40*cm2; 600 theVec->PutValue(k,logene,G4Log(aValue)); 601 } 602 } 603 604 if (fVerboseLevel > 2) 605 { 606 G4cout << "G4PenelopePhotoElectricModel: 607 << Z << G4endl; 608 } 609 610 fLogAtomicShellXS[Z] = thePhysicsTable; 611 612 file.close(); 613 return; 614 } 615 616 //....oooOO0OOooo........oooOO0OOooo........oo 617 618 std::size_t G4PenelopePhotoElectricModel::GetN 619 { 620 if (!IsMaster()) 621 //Should not be here! 622 G4Exception("G4PenelopePhotoElectricModel: 623 "em0100",FatalException,"Worker thread in 624 625 //read data files 626 if (!fLogAtomicShellXS[Z]) 627 ReadDataFile(Z); 628 //now it should be ok 629 if (!fLogAtomicShellXS[Z]) 630 { 631 G4ExceptionDescription ed; 632 ed << "Cannot find shell cross section 633 G4Exception("G4PenelopePhotoElectricMod 634 "em2038",FatalException,ed); 635 } 636 //one vector is allocated for the _total_ cr 637 std::size_t nEntries = fLogAtomicShellXS[Z]- 638 return (nEntries-1); 639 } 640 641 //....oooOO0OOooo........oooOO0OOooo........oo 642 643 G4double G4PenelopePhotoElectricModel::GetShel 644 { 645 //this forces also the loading of the data 646 std::size_t entries = GetNumberOfShellXS(Z); 647 648 if (shellID >= entries) 649 { 650 G4cout << "Element Z=" << Z << " has dat 651 G4cout << "so shellID should be from 0 t 652 return 0; 653 } 654 655 G4PhysicsTable* theTable = fLogAtomicShellX 656 //[0] is the total XS, shellID is in the ele 657 G4PhysicsFreeVector* totalXSLog = (G4Physics 658 659 if (!totalXSLog) 660 { 661 G4Exception("G4PenelopePhotoElectricMod 662 "em2039",FatalException, 663 "Unable to retrieve the total cross sec 664 return 0; 665 } 666 G4double logene = G4Log(energy); 667 G4double logXS = totalXSLog->Value(logene); 668 G4double cross = G4Exp(logXS); 669 if (cross < 2e-40*cm2) cross = 0; 670 return cross; 671 } 672 673 //....oooOO0OOooo........oooOO0OOooo........oo 674 675 G4String G4PenelopePhotoElectricModel::WriteTa 676 { 677 G4String theShell = "outer shell"; 678 if (shellID == 0) 679 theShell = "K"; 680 else if (shellID == 1) 681 theShell = "L1"; 682 else if (shellID == 2) 683 theShell = "L2"; 684 else if (shellID == 3) 685 theShell = "L3"; 686 else if (shellID == 4) 687 theShell = "M1"; 688 else if (shellID == 5) 689 theShell = "M2"; 690 else if (shellID == 6) 691 theShell = "M3"; 692 else if (shellID == 7) 693 theShell = "M4"; 694 else if (shellID == 8) 695 theShell = "M5"; 696 697 return theShell; 698 } 699 700 //....oooOO0OOooo........oooOO0OOooo........oo 701 702 void G4PenelopePhotoElectricModel::SetParticle 703 { 704 if(!fParticle) { 705 fParticle = p; 706 } 707 } 708 709 //....oooOO0OOooo........oooOO0OOooo........oo 710 711 std::size_t G4PenelopePhotoElectricModel::Sele 712 { 713 G4double logEnergy = G4Log(energy); 714 715 //Check if data have been read (it should be 716 if (!fLogAtomicShellXS[Z]) 717 { 718 G4ExceptionDescription ed; 719 ed << "Cannot find shell cross section 720 G4Exception("G4PenelopePhotoElectricMod 721 "em2038",FatalException,ed); 722 } 723 724 G4PhysicsTable* theTable = fLogAtomicShellX 725 726 G4double sum = 0; 727 G4PhysicsFreeVector* totalXSLog = (G4Physics 728 G4double logXS = totalXSLog->Value(logEnergy 729 G4double totalXS = G4Exp(logXS); 730 731 //Notice: totalXS is the total cross section 732 //the sum of partialXS's, since these includ 733 // 734 // Therefore, here one have to consider the 735 // an outer shell. Conventionally, it is ind 736 // 737 G4double random = G4UniformRand()*totalXS; 738 739 for (std::size_t k=1;k<theTable->entries();+ 740 { 741 //Add one shell 742 G4PhysicsFreeVector* partialXSLog = (G4P 743 G4double logXSLocal = partialXSLog->Valu 744 G4double partialXS = G4Exp(logXSLocal); 745 sum += partialXS; 746 if (random <= sum) 747 return k-1; 748 } 749 //none of the shells K, L, M: return outer s 750 return 9; 751 } 752