00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include "G4PenelopePhotoElectricModel.hh"
00042 #include "G4PhysicalConstants.hh"
00043 #include "G4SystemOfUnits.hh"
00044 #include "G4ParticleDefinition.hh"
00045 #include "G4MaterialCutsCouple.hh"
00046 #include "G4DynamicParticle.hh"
00047 #include "G4PhysicsTable.hh"
00048 #include "G4PhysicsFreeVector.hh"
00049 #include "G4ElementTable.hh"
00050 #include "G4Element.hh"
00051 #include "G4AtomicTransitionManager.hh"
00052 #include "G4AtomicShell.hh"
00053 #include "G4Gamma.hh"
00054 #include "G4Electron.hh"
00055 #include "G4LossTableManager.hh"
00056
00057
00058
00059
00060 G4PenelopePhotoElectricModel::G4PenelopePhotoElectricModel(const G4ParticleDefinition*,
00061 const G4String& nam)
00062 :G4VEmModel(nam),fParticleChange(0),isInitialised(false),fAtomDeexcitation(0),
00063 logAtomicShellXS(0)
00064 {
00065 fIntrinsicLowEnergyLimit = 100.0*eV;
00066 fIntrinsicHighEnergyLimit = 100.0*GeV;
00067
00068 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
00069
00070 verboseLevel= 0;
00071
00072
00073
00074
00075
00076
00077
00078
00079 SetDeexcitationFlag(true);
00080
00081 fTransitionManager = G4AtomicTransitionManager::Instance();
00082 }
00083
00084
00085
00086 G4PenelopePhotoElectricModel::~G4PenelopePhotoElectricModel()
00087 {
00088 std::map <G4int,G4PhysicsTable*>::iterator i;
00089 if (logAtomicShellXS)
00090 {
00091 for (i=logAtomicShellXS->begin();i != logAtomicShellXS->end();i++)
00092 {
00093 G4PhysicsTable* tab = i->second;
00094 tab->clearAndDestroy();
00095 delete tab;
00096 }
00097 }
00098 delete logAtomicShellXS;
00099 }
00100
00101
00102
00103 void G4PenelopePhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
00104 const G4DataVector& cuts)
00105 {
00106 if (verboseLevel > 3)
00107 G4cout << "Calling G4PenelopePhotoElectricModel::Initialise()" << G4endl;
00108
00109
00110 if (!logAtomicShellXS)
00111 logAtomicShellXS = new std::map<G4int,G4PhysicsTable*>;
00112
00113 InitialiseElementSelectors(particle,cuts);
00114
00115 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00116
00117 if (!fAtomDeexcitation)
00118 {
00119 G4cout << G4endl;
00120 G4cout << "WARNING from G4PenelopePhotoElectricModel " << G4endl;
00121 G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
00122 G4cout << "any fluorescence/Auger emission." << G4endl;
00123 G4cout << "Please make sure this is intended" << G4endl;
00124 }
00125
00126 if (verboseLevel > 0) {
00127 G4cout << "Penelope Photo-Electric model v2008 is initialized " << G4endl
00128 << "Energy range: "
00129 << LowEnergyLimit() / MeV << " MeV - "
00130 << HighEnergyLimit() / GeV << " GeV";
00131 }
00132
00133 if(isInitialised) return;
00134 fParticleChange = GetParticleChangeForGamma();
00135 isInitialised = true;
00136
00137 }
00138
00139
00140
00141 G4double G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom(
00142 const G4ParticleDefinition*,
00143 G4double energy,
00144 G4double Z, G4double,
00145 G4double, G4double)
00146 {
00147
00148
00149
00150
00151 if (verboseLevel > 3)
00152 G4cout << "Calling ComputeCrossSectionPerAtom() of G4PenelopePhotoElectricModel" << G4endl;
00153
00154 G4int iZ = (G4int) Z;
00155
00156
00157 if (!logAtomicShellXS->count(iZ))
00158 ReadDataFile(iZ);
00159
00160 if (!logAtomicShellXS->count(iZ))
00161 G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
00162 "em2038",FatalException,
00163 "Unable to retrieve the shell cross section table");
00164
00165 G4double cross = 0;
00166
00167 G4PhysicsTable* theTable = logAtomicShellXS->find(iZ)->second;
00168 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
00169
00170 if (!totalXSLog)
00171 {
00172 G4Exception("G4PenelopePhotoElectricModel::ComputeCrossSectionPerAtom()",
00173 "em2039",FatalException,
00174 "Unable to retrieve the total cross section table");
00175 return 0;
00176 }
00177 G4double logene = std::log(energy);
00178 G4double logXS = totalXSLog->Value(logene);
00179 cross = std::exp(logXS);
00180
00181 if (verboseLevel > 2)
00182 G4cout << "Photoelectric cross section at " << energy/MeV << " MeV for Z=" << Z <<
00183 " = " << cross/barn << " barn" << G4endl;
00184 return cross;
00185 }
00186
00187
00188
00189 void G4PenelopePhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00190 const G4MaterialCutsCouple* couple,
00191 const G4DynamicParticle* aDynamicGamma,
00192 G4double,
00193 G4double)
00194 {
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209 if (verboseLevel > 3)
00210 G4cout << "Calling SamplingSecondaries() of G4PenelopePhotoElectricModel" << G4endl;
00211
00212 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00213
00214
00215 fParticleChange->ProposeTrackStatus(fStopAndKill);
00216 fParticleChange->SetProposedKineticEnergy(0.);
00217
00218 if (photonEnergy <= fIntrinsicLowEnergyLimit)
00219 {
00220 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00221 return ;
00222 }
00223
00224 G4ParticleMomentum photonDirection = aDynamicGamma->GetMomentumDirection();
00225
00226
00227 if (verboseLevel > 2)
00228 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
00229
00230
00231 const G4Element* anElement =
00232 SelectRandomAtom(couple,G4Gamma::GammaDefinition(),photonEnergy);
00233 G4int Z = (G4int) anElement->GetZ();
00234 if (verboseLevel > 2)
00235 G4cout << "Selected " << anElement->GetName() << G4endl;
00236
00237
00238
00239
00240
00241
00242
00243 size_t shellIndex = SelectRandomShell(Z,photonEnergy);
00244
00245 if (verboseLevel > 2)
00246 G4cout << "Selected shell " << shellIndex << " of element " << anElement->GetName() << G4endl;
00247
00248
00249 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
00250
00251
00252
00253
00254
00255
00256
00257 size_t numberOfShells = (size_t) transitionManager->NumberOfShells(Z);
00258 if (shellIndex >= numberOfShells)
00259 shellIndex = numberOfShells-1;
00260
00261 const G4AtomicShell* shell = fTransitionManager->Shell(Z,shellIndex);
00262 G4double bindingEnergy = shell->BindingEnergy();
00263
00264
00265
00266
00267
00268
00269
00270 if (shellIndex == 9)
00271 bindingEnergy = 0.*eV;
00272
00273
00274 G4double localEnergyDeposit = 0.0;
00275 G4double cosTheta = 1.0;
00276
00277
00278 G4double eKineticEnergy = photonEnergy - bindingEnergy;
00279
00280
00281
00282 if (eKineticEnergy > 0.)
00283 {
00284
00285
00286 cosTheta = SampleElectronDirection(eKineticEnergy);
00287 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
00288 G4double phi = twopi * G4UniformRand() ;
00289 G4double dirx = sinTheta * std::cos(phi);
00290 G4double diry = sinTheta * std::sin(phi);
00291 G4double dirz = cosTheta ;
00292 G4ThreeVector electronDirection(dirx,diry,dirz);
00293 electronDirection.rotateUz(photonDirection);
00294 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),
00295 electronDirection,
00296 eKineticEnergy);
00297 fvect->push_back(electron);
00298 }
00299 else
00300 bindingEnergy = photonEnergy;
00301
00302
00303 G4double energyInFluorescence = 0;
00304 G4double energyInAuger = 0;
00305
00306
00307
00308
00309 if (fAtomDeexcitation && shellIndex<9)
00310 {
00311 G4int index = couple->GetIndex();
00312 if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
00313 {
00314 size_t nBefore = fvect->size();
00315 fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
00316 size_t nAfter = fvect->size();
00317
00318 if (nAfter > nBefore)
00319 {
00320 for (size_t j=nBefore;j<nAfter;j++)
00321 {
00322 G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
00323 bindingEnergy -= itsEnergy;
00324 if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
00325 energyInFluorescence += itsEnergy;
00326 else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
00327 energyInAuger += itsEnergy;
00328
00329 }
00330 }
00331
00332 }
00333 }
00334
00335
00336 localEnergyDeposit += bindingEnergy;
00337
00338 if (localEnergyDeposit < 0)
00339 {
00340 G4cout << "WARNING - "
00341 << "G4PenelopePhotoElectricModel::SampleSecondaries() - Negative energy deposit"
00342 << G4endl;
00343 localEnergyDeposit = 0;
00344 }
00345
00346 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
00347
00348 if (verboseLevel > 1)
00349 {
00350 G4cout << "-----------------------------------------------------------" << G4endl;
00351 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
00352 G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " <<
00353 anElement->GetName() << G4endl;
00354 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
00355 G4cout << "-----------------------------------------------------------" << G4endl;
00356 if (eKineticEnergy)
00357 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
00358 if (energyInFluorescence)
00359 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
00360 if (energyInAuger)
00361 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
00362 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
00363 G4cout << "Total final state: " <<
00364 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV <<
00365 " keV" << G4endl;
00366 G4cout << "-----------------------------------------------------------" << G4endl;
00367 }
00368 if (verboseLevel > 0)
00369 {
00370 G4double energyDiff =
00371 std::fabs(eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger-photonEnergy);
00372 if (energyDiff > 0.05*keV)
00373 {
00374 G4cout << "Warning from G4PenelopePhotoElectric: problem with energy conservation: " <<
00375 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV
00376 << " keV (final) vs. " <<
00377 photonEnergy/keV << " keV (initial)" << G4endl;
00378 G4cout << "-----------------------------------------------------------" << G4endl;
00379 G4cout << "Energy balance from G4PenelopePhotoElectric" << G4endl;
00380 G4cout << "Selected shell: " << WriteTargetShell(shellIndex) << " of element " <<
00381 anElement->GetName() << G4endl;
00382 G4cout << "Incoming photon energy: " << photonEnergy/keV << " keV" << G4endl;
00383 G4cout << "-----------------------------------------------------------" << G4endl;
00384 if (eKineticEnergy)
00385 G4cout << "Outgoing electron " << eKineticEnergy/keV << " keV" << G4endl;
00386 if (energyInFluorescence)
00387 G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
00388 if (energyInAuger)
00389 G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
00390 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
00391 G4cout << "Total final state: " <<
00392 (eKineticEnergy+energyInFluorescence+localEnergyDeposit+energyInAuger)/keV <<
00393 " keV" << G4endl;
00394 G4cout << "-----------------------------------------------------------" << G4endl;
00395 }
00396 }
00397 }
00398
00399
00400
00401 G4double G4PenelopePhotoElectricModel::SampleElectronDirection(G4double energy)
00402 {
00403 G4double costheta = 1.0;
00404 if (energy>1*GeV) return costheta;
00405
00406
00407
00408
00409 G4double gamma = 1.0 + energy/electron_mass_c2;
00410 G4double gamma2 = gamma*gamma;
00411 G4double beta = std::sqrt((gamma2-1.0)/gamma2);
00412
00413
00414
00415 G4double ac = (1.0/beta) - 1.0;
00416 G4double a1 = 0.5*beta*gamma*(gamma-1.0)*(gamma-2.0);
00417 G4double a2 = ac + 2.0;
00418 G4double gtmax = 2.0*(a1 + 1.0/ac);
00419
00420 G4double tsam = 0;
00421 G4double gtr = 0;
00422
00423
00424
00425
00426 do{
00427 G4double rand = G4UniformRand();
00428 tsam = 2.0*ac * (2.0*rand + a2*std::sqrt(rand)) / (a2*a2 - 4.0*rand);
00429 gtr = (2.0 - tsam) * (a1 + 1.0/(ac+tsam));
00430 }while(G4UniformRand()*gtmax > gtr);
00431 costheta = 1.0-tsam;
00432
00433
00434 return costheta;
00435 }
00436
00437
00438
00439 void G4PenelopePhotoElectricModel::ReadDataFile(G4int Z)
00440 {
00441 if (verboseLevel > 2)
00442 {
00443 G4cout << "G4PenelopePhotoElectricModel::ReadDataFile()" << G4endl;
00444 G4cout << "Going to read PhotoElectric data files for Z=" << Z << G4endl;
00445 }
00446
00447 char* path = getenv("G4LEDATA");
00448 if (!path)
00449 {
00450 G4String excep = "G4PenelopePhotoElectricModel - G4LEDATA environment variable not set!";
00451 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
00452 "em0006",FatalException,excep);
00453 return;
00454 }
00455
00456
00457
00458
00459 std::ostringstream ost;
00460 if (Z>9)
00461 ost << path << "/penelope/photoelectric/pdgph" << Z << ".p08";
00462 else
00463 ost << path << "/penelope/photoelectric/pdgph0" << Z << ".p08";
00464 std::ifstream file(ost.str().c_str());
00465 if (!file.is_open())
00466 {
00467 G4String excep = "G4PenelopePhotoElectricModel - data file " + G4String(ost.str()) + " not found!";
00468 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
00469 "em0003",FatalException,excep);
00470 }
00471
00472
00473 size_t ndata=0;
00474 G4String line;
00475 while( getline(file, line) )
00476 ndata++;
00477 ndata -= 1;
00478
00479
00480 file.clear();
00481 file.close();
00482 file.open(ost.str().c_str());
00483
00484 G4int readZ =0;
00485 size_t nShells= 0;
00486 file >> readZ >> nShells;
00487
00488 if (verboseLevel > 3)
00489 G4cout << "Element Z=" << Z << " , nShells = " << nShells << G4endl;
00490
00491
00492 if (readZ != Z || nShells <= 0 || nShells > 50)
00493 {
00494 G4ExceptionDescription ed;
00495 ed << "Corrupted data file for Z=" << Z << G4endl;
00496 G4Exception("G4PenelopePhotoElectricModel::ReadDataFile()",
00497 "em0005",FatalException,ed);
00498 return;
00499 }
00500 G4PhysicsTable* thePhysicsTable = new G4PhysicsTable();
00501
00502
00503
00504
00505
00506
00507
00508 for (size_t i=0;i<nShells+1;i++)
00509 thePhysicsTable->push_back(new G4PhysicsFreeVector(ndata));
00510
00511 size_t k =0;
00512 for (k=0;k<ndata && !file.eof();k++)
00513 {
00514 G4double energy = 0;
00515 G4double aValue = 0;
00516 file >> energy ;
00517 energy *= eV;
00518 G4double logene = std::log(energy);
00519
00520 for (size_t i=0;i<nShells+1;i++)
00521 {
00522 file >> aValue;
00523 aValue *= barn;
00524 G4PhysicsFreeVector* theVec = (G4PhysicsFreeVector*) ((*thePhysicsTable)[i]);
00525 if (aValue < 1e-40*cm2)
00526 aValue = 1e-40*cm2;
00527 theVec->PutValue(k,logene,std::log(aValue));
00528 }
00529 }
00530
00531 if (verboseLevel > 2)
00532 {
00533 G4cout << "G4PenelopePhotoElectricModel: read " << k << " points for element Z = "
00534 << Z << G4endl;
00535 }
00536
00537 logAtomicShellXS->insert(std::make_pair(Z,thePhysicsTable));
00538
00539 file.close();
00540 return;
00541 }
00542
00543
00544
00545 size_t G4PenelopePhotoElectricModel::SelectRandomShell(G4int Z,G4double energy)
00546 {
00547 G4double logEnergy = std::log(energy);
00548
00549
00550 if (!logAtomicShellXS->count(Z))
00551 {
00552 G4ExceptionDescription ed;
00553 ed << "Cannot find shell cross section data for Z=" << Z << G4endl;
00554 G4Exception("G4PenelopePhotoElectricModel::SelectRandomShell()",
00555 "em2038",FatalException,ed);
00556 }
00557
00558 size_t shellIndex = 0;
00559
00560 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
00561
00562 G4DataVector* tempVector = new G4DataVector();
00563
00564 G4double sum = 0;
00565
00566
00567 tempVector->push_back(sum);
00568
00569 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[0];
00570 G4double logXS = totalXSLog->Value(logEnergy);
00571 G4double totalXS = std::exp(logXS);
00572
00573
00574
00575
00576
00577
00578
00579
00580 for (size_t k=1;k<theTable->entries();k++)
00581 {
00582 G4PhysicsFreeVector* partialXSLog = (G4PhysicsFreeVector*) (*theTable)[k];
00583 G4double logXSLocal = partialXSLog->Value(logEnergy);
00584 G4double partialXS = std::exp(logXSLocal);
00585 sum += partialXS;
00586 tempVector->push_back(sum);
00587 }
00588
00589 tempVector->push_back(totalXS);
00590
00591 G4double random = G4UniformRand()*totalXS;
00592
00593
00594
00595
00596
00597
00598
00599
00600
00601 size_t lowerBound = 0;
00602 size_t upperBound = tempVector->size()-1;
00603 while (lowerBound <= upperBound)
00604 {
00605 size_t midBin = (lowerBound + upperBound)/2;
00606 if( random < (*tempVector)[midBin])
00607 upperBound = midBin-1;
00608 else
00609 lowerBound = midBin+1;
00610 }
00611
00612 shellIndex = upperBound;
00613
00614 delete tempVector;
00615 return shellIndex;
00616 }
00617
00618
00619
00620 size_t G4PenelopePhotoElectricModel::GetNumberOfShellXS(G4int Z)
00621 {
00622
00623 if (!logAtomicShellXS->count(Z))
00624 ReadDataFile(Z);
00625
00626 if (!logAtomicShellXS->count(Z))
00627 {
00628 G4ExceptionDescription ed;
00629 ed << "Cannot find shell cross section data for Z=" << Z << G4endl;
00630 G4Exception("G4PenelopePhotoElectricModel::GetNumberOfShellXS()",
00631 "em2038",FatalException,ed);
00632 }
00633
00634 size_t nEntries = logAtomicShellXS->find(Z)->second->entries();
00635 return (nEntries-1);
00636 }
00637
00638
00639
00640 G4double G4PenelopePhotoElectricModel::GetShellCrossSection(G4int Z,size_t shellID,G4double energy)
00641 {
00642
00643 size_t entries = GetNumberOfShellXS(Z);
00644
00645 if (shellID >= entries)
00646 {
00647 G4cout << "Element Z=" << Z << " has data for " << entries << " shells only" << G4endl;
00648 G4cout << "so shellID should be from 0 to " << entries-1 << G4endl;
00649 return 0;
00650 }
00651
00652 G4PhysicsTable* theTable = logAtomicShellXS->find(Z)->second;
00653
00654 G4PhysicsFreeVector* totalXSLog = (G4PhysicsFreeVector*) (*theTable)[shellID+1];
00655
00656 if (!totalXSLog)
00657 {
00658 G4Exception("G4PenelopePhotoElectricModel::GetShellCrossSection()",
00659 "em2039",FatalException,
00660 "Unable to retrieve the total cross section table");
00661 return 0;
00662 }
00663 G4double logene = std::log(energy);
00664 G4double logXS = totalXSLog->Value(logene);
00665 G4double cross = std::exp(logXS);
00666 if (cross < 2e-40*cm2) cross = 0;
00667 return cross;
00668 }
00669
00670
00671
00672 G4String G4PenelopePhotoElectricModel::WriteTargetShell(size_t shellID)
00673 {
00674 G4String theShell = "outer shell";
00675 if (shellID == 0)
00676 theShell = "K";
00677 else if (shellID == 1)
00678 theShell = "L1";
00679 else if (shellID == 2)
00680 theShell = "L2";
00681 else if (shellID == 3)
00682 theShell = "L3";
00683 else if (shellID == 4)
00684 theShell = "M1";
00685 else if (shellID == 5)
00686 theShell = "M2";
00687 else if (shellID == 6)
00688 theShell = "M3";
00689 else if (shellID == 7)
00690 theShell = "M4";
00691 else if (shellID == 8)
00692 theShell = "M5";
00693
00694 return theShell;
00695 }