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
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 #include "G4UAtomicDeexcitation.hh"
00053 #include "G4PhysicalConstants.hh"
00054 #include "G4SystemOfUnits.hh"
00055 #include "Randomize.hh"
00056 #include "G4Gamma.hh"
00057 #include "G4AtomicTransitionManager.hh"
00058 #include "G4FluoTransition.hh"
00059 #include "G4Electron.hh"
00060 #include "G4Positron.hh"
00061 #include "G4Proton.hh"
00062 #include "G4Alpha.hh"
00063
00064 #include "G4teoCrossSection.hh"
00065 #include "G4empCrossSection.hh"
00066 #include "G4PenelopeIonisationCrossSection.hh"
00067 #include "G4LivermoreIonisationCrossSection.hh"
00068 #include "G4EmCorrections.hh"
00069 #include "G4LossTableManager.hh"
00070 #include "G4Material.hh"
00071 #include "G4AtomicShells.hh"
00072
00073 using namespace std;
00074
00075 G4UAtomicDeexcitation::G4UAtomicDeexcitation():
00076 G4VAtomDeexcitation("UAtomDeexcitation"),
00077 minGammaEnergy(DBL_MAX),
00078 minElectronEnergy(DBL_MAX),
00079 emcorr(0)
00080 {
00081 PIXEshellCS = 0;
00082 ePIXEshellCS = 0;
00083 emcorr = G4LossTableManager::Instance()->EmCorrections();
00084 theElectron = G4Electron::Electron();
00085 thePositron = G4Positron::Positron();
00086 transitionManager = 0;
00087 anaPIXEshellCS = 0;
00088 }
00089
00090 G4UAtomicDeexcitation::~G4UAtomicDeexcitation()
00091 {
00092 delete PIXEshellCS;
00093 delete anaPIXEshellCS;
00094 delete ePIXEshellCS;
00095 }
00096
00097 void G4UAtomicDeexcitation::InitialiseForNewRun()
00098 {
00099 if(!IsFluoActive()) { return; }
00100 transitionManager = G4AtomicTransitionManager::Instance();
00101 if(IsPIXEActive()) {
00102 G4cout << G4endl;
00103 G4cout << "### === G4UAtomicDeexcitation::InitialiseForNewRun()" << G4endl;
00104 anaPIXEshellCS = new G4teoCrossSection("Analytical");
00105
00106 }
00107 else {return;}
00108
00109
00110 if (PIXECrossSectionModel() == "" ||
00111 PIXECrossSectionModel() == "Empirical" ||
00112 PIXECrossSectionModel() == "empirical")
00113 {
00114 SetPIXECrossSectionModel("Empirical");
00115 }
00116 else if (PIXECrossSectionModel() == "ECPSSR_Analytical" ||
00117 PIXECrossSectionModel() == "Analytical" ||
00118 PIXECrossSectionModel() == "analytical")
00119 {
00120 SetPIXECrossSectionModel("Analytical");
00121 }
00122 else if (PIXECrossSectionModel() == "ECPSSR_FormFactor" ||
00123 PIXECrossSectionModel() == "ECPSSR_Tabulated" ||
00124 PIXECrossSectionModel() == "Analytical_Tabulated")
00125 {
00126 SetPIXECrossSectionModel("ECPSSR_FormFactor");
00127 }
00128 else
00129 {
00130 G4cout << "### G4UAtomicDeexcitation::InitialiseForNewRun WARNING "
00131 << G4endl;
00132 G4cout << " PIXE cross section name " << PIXECrossSectionModel()
00133 << " is unknown, Analytical cross section will be used" << G4endl;
00134 SetPIXECrossSectionModel("Analytical");
00135 }
00136
00137
00138 if(PIXEshellCS)
00139 {
00140 if(PIXECrossSectionModel() != PIXEshellCS->GetName())
00141 {
00142 delete PIXEshellCS;
00143 PIXEshellCS = 0;
00144 }
00145 }
00146
00147
00148 if(!PIXEshellCS) {
00149 if (PIXECrossSectionModel() == "Empirical")
00150 {
00151 PIXEshellCS = new G4empCrossSection("Empirical");
00152 }
00153
00154 if (PIXECrossSectionModel() == "ECPSSR_FormFactor")
00155 {
00156 PIXEshellCS = new G4teoCrossSection("ECPSSR_FormFactor");
00157 }
00158 }
00159
00160
00161
00162
00163 if (PIXEElectronCrossSectionModel() == "" ||
00164 PIXEElectronCrossSectionModel() == "Livermore")
00165 {
00166 SetPIXEElectronCrossSectionModel("Livermore");
00167 }
00168 else if (PIXEElectronCrossSectionModel() == "ProtonAnalytical" ||
00169 PIXEElectronCrossSectionModel() == "Analytical" ||
00170 PIXEElectronCrossSectionModel() == "analytical")
00171 {
00172 SetPIXEElectronCrossSectionModel("ProtonAnalytical");
00173 }
00174 else if (PIXEElectronCrossSectionModel() == "ProtonEmpirical" ||
00175 PIXEElectronCrossSectionModel() == "Empirical" ||
00176 PIXEElectronCrossSectionModel() == "empirical")
00177 {
00178 SetPIXEElectronCrossSectionModel("ProtonEmpirical");
00179 }
00180 else if (PIXEElectronCrossSectionModel() == "Penelope")
00181 SetPIXEElectronCrossSectionModel("Penelope");
00182 else
00183 {
00184 G4cout << "### G4UAtomicDeexcitation::InitialiseForNewRun WARNING "
00185 << G4endl;
00186 G4cout << " PIXE e- cross section name " << PIXEElectronCrossSectionModel()
00187 << " is unknown, PIXE is disabled" << G4endl;
00188 SetPIXEElectronCrossSectionModel("Livermore");
00189 }
00190
00191
00192 if(ePIXEshellCS)
00193 {
00194 if(PIXEElectronCrossSectionModel() != ePIXEshellCS->GetName())
00195 {
00196 delete ePIXEshellCS;
00197 ePIXEshellCS = 0;
00198 }
00199 }
00200
00201
00202 if(!ePIXEshellCS)
00203 {
00204 if(PIXEElectronCrossSectionModel() == "Empirical")
00205 {
00206 ePIXEshellCS = new G4empCrossSection("Empirical");
00207 }
00208
00209 else if(PIXEElectronCrossSectionModel() == "Analytical")
00210 {
00211 ePIXEshellCS = new G4teoCrossSection("Analytical");
00212 }
00213
00214 else if(PIXEElectronCrossSectionModel() == "Livermore")
00215 {
00216 ePIXEshellCS = new G4LivermoreIonisationCrossSection();
00217 }
00218 else if (PIXEElectronCrossSectionModel() == "Penelope")
00219 {
00220 ePIXEshellCS = new G4PenelopeIonisationCrossSection();
00221 }
00222 }
00223 }
00224
00225 void G4UAtomicDeexcitation::InitialiseForExtraAtom(G4int )
00226 {}
00227
00228 const G4AtomicShell* G4UAtomicDeexcitation::GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)
00229 {
00230 return transitionManager->Shell(Z, size_t(shell));
00231 }
00232
00233 void G4UAtomicDeexcitation::GenerateParticles(
00234 std::vector<G4DynamicParticle*>* vectorOfParticles,
00235 const G4AtomicShell* atomicShell,
00236 G4int Z,
00237 G4double gammaCut,
00238 G4double eCut)
00239 {
00240
00241
00242 G4int givenShellId = atomicShell->ShellId();
00243
00244 minGammaEnergy = gammaCut;
00245 minElectronEnergy = eCut;
00246
00247
00248
00249
00250
00251 G4DynamicParticle* aParticle=0;
00252 G4int provShellId = 0;
00253 G4int counter = 0;
00254
00255
00256
00257 if (Z>5 && Z<100) {
00258
00259
00260
00261 do
00262 {
00263 if (counter == 0)
00264
00265
00266 {
00267 provShellId = SelectTypeOfTransition(Z, givenShellId);
00268
00269 if ( provShellId >0)
00270 {
00271 aParticle = GenerateFluorescence(Z,givenShellId,provShellId);
00272
00273 }
00274 else if ( provShellId == -1)
00275 {
00276
00277 aParticle = GenerateAuger(Z, givenShellId);
00278
00279 }
00280 else
00281 {
00282 G4Exception("G4UAtomicDeexcitation::GenerateParticles()","de0002",JustWarning, "Energy deposited locally");
00283 }
00284 }
00285 else
00286
00287
00288 {
00289 provShellId = SelectTypeOfTransition(Z,newShellId);
00290 if (provShellId >0)
00291 {
00292 aParticle = GenerateFluorescence(Z,newShellId,provShellId);
00293
00294 }
00295 else if ( provShellId == -1)
00296 {
00297
00298 aParticle = GenerateAuger(Z, newShellId);
00299
00300 }
00301 else
00302 {
00303 G4Exception("G4UAtomicDeexcitation::GenerateParticles()","de0002",JustWarning, "Energy deposited locally");
00304 }
00305 }
00306 counter++;
00307 if (aParticle != 0)
00308 {
00309 vectorOfParticles->push_back(aParticle);
00310
00311 }
00312 else {provShellId = -2;}
00313 }
00314 while (provShellId > -2);
00315 }
00316 else
00317 {
00318 G4Exception("G4UAtomicDeexcitation::GenerateParticles()","de0001",JustWarning, "Energy deposited locally");
00319 }
00320
00321
00322
00323 }
00324
00325 G4double
00326 G4UAtomicDeexcitation::GetShellIonisationCrossSectionPerAtom(
00327 const G4ParticleDefinition* pdef,
00328 G4int Z,
00329 G4AtomicShellEnumerator shellEnum,
00330 G4double kineticEnergy,
00331 const G4Material* mat)
00332 {
00333
00334
00335
00336
00337
00338 G4double xsec = 0.0;
00339 if(Z > 93 || Z < 6 ) { return xsec; }
00340 G4int idx = G4int(shellEnum);
00341 if(idx >= G4AtomicShells::GetNumberOfShells(Z)) { return xsec; }
00342
00343
00344 if(pdef == theElectron || pdef == thePositron) {
00345 xsec = ePIXEshellCS->CrossSection(Z,shellEnum,kineticEnergy,0.0,mat);
00346 return xsec;
00347 }
00348
00349 G4double mass = pdef->GetPDGMass();
00350 G4double escaled = kineticEnergy;
00351 G4double q2 = 0.0;
00352
00353
00354 if ((pdef->GetParticleName() != "proton" && pdef->GetParticleName() != "alpha" ) )
00355 {
00356 mass = proton_mass_c2;
00357 escaled = kineticEnergy*mass/(pdef->GetPDGMass());
00358
00359 if(mat) {
00360 q2 = emcorr->EffectiveChargeSquareRatio(pdef,mat,kineticEnergy);
00361 } else {
00362 G4double q = pdef->GetPDGCharge()/eplus;
00363 q2 = q*q;
00364 }
00365 }
00366
00367 if(PIXEshellCS) { xsec = PIXEshellCS->CrossSection(Z,shellEnum,escaled,mass,mat); }
00368 if(xsec < 1e-100) {
00369
00370 xsec = anaPIXEshellCS->CrossSection(Z,shellEnum,escaled,mass,mat);
00371
00372 }
00373
00374 if (q2) {xsec *= q2;}
00375
00376 return xsec;
00377 }
00378
00379 void G4UAtomicDeexcitation::SetCutForSecondaryPhotons(G4double cut)
00380 {
00381 minGammaEnergy = cut;
00382 }
00383
00384 void G4UAtomicDeexcitation::SetCutForAugerElectrons(G4double cut)
00385 {
00386 minElectronEnergy = cut;
00387 }
00388
00389 G4double
00390 G4UAtomicDeexcitation::ComputeShellIonisationCrossSectionPerAtom(
00391 const G4ParticleDefinition* p,
00392 G4int Z,
00393 G4AtomicShellEnumerator shell,
00394 G4double kinE,
00395 const G4Material* mat)
00396 {
00397 return GetShellIonisationCrossSectionPerAtom(p,Z,shell,kinE,mat);
00398 }
00399
00400 G4int G4UAtomicDeexcitation::SelectTypeOfTransition(G4int Z, G4int shellId)
00401 {
00402 if (shellId <=0 ) {
00403 G4Exception("G4UAtomicDeexcitation::SelecttypeOfTransition()","de0002",JustWarning, "Energy deposited locally");
00404 return 0;
00405 }
00406
00407
00408 G4int provShellId = -1;
00409 G4int shellNum = 0;
00410 G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
00411
00412 const G4FluoTransition* refShell = transitionManager->ReachableShell(Z,maxNumOfShells-1);
00413
00414
00415
00416
00417 if ( shellId <= refShell->FinalShellId())
00418 {
00419 while (shellId != transitionManager->ReachableShell(Z,shellNum)->FinalShellId())
00420 {
00421 if(shellNum ==maxNumOfShells-1)
00422 {
00423 break;
00424 }
00425 shellNum++;
00426 }
00427 G4int transProb = 0;
00428
00429 G4double partialProb = G4UniformRand();
00430 G4double partSum = 0;
00431 const G4FluoTransition* aShell = transitionManager->ReachableShell(Z,shellNum);
00432 G4int trSize = (aShell->TransitionProbabilities()).size();
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442 while(transProb < trSize){
00443
00444 partSum += aShell->TransitionProbability(transProb);
00445
00446 if(partialProb <= partSum)
00447 {
00448 provShellId = aShell->OriginatingShellId(transProb);
00449
00450
00451 break;
00452 }
00453 transProb++;
00454 }
00455
00456
00457
00458 }
00459 else
00460 {
00461 provShellId = -1;
00462 }
00463
00464
00465 return provShellId;
00466 }
00467
00468 G4DynamicParticle*
00469 G4UAtomicDeexcitation::GenerateFluorescence(G4int Z, G4int shellId,
00470 G4int provShellId )
00471 {
00472
00473
00474
00475 if (shellId <=0 )
00476
00477 {
00478 G4Exception("G4UAtomicDeexcitation::GenerateFluorescence()","de0002",JustWarning, "Energy deposited locally");
00479 return 0;
00480 }
00481
00482
00483
00484 G4double newcosTh = 1.-2.*G4UniformRand();
00485 G4double newsinTh = std::sqrt((1.-newcosTh)*(1. + newcosTh));
00486 G4double newPhi = twopi*G4UniformRand();
00487
00488 G4double xDir = newsinTh*std::sin(newPhi);
00489 G4double yDir = newsinTh*std::cos(newPhi);
00490 G4double zDir = newcosTh;
00491
00492 G4ThreeVector newGammaDirection(xDir,yDir,zDir);
00493
00494 G4int shellNum = 0;
00495 G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
00496
00497
00498 while (shellId != transitionManager->
00499 ReachableShell(Z,shellNum)->FinalShellId())
00500 {
00501 if(shellNum == maxNumOfShells-1)
00502 {
00503 break;
00504 }
00505 shellNum++;
00506 }
00507
00508 size_t transitionSize = transitionManager->
00509 ReachableShell(Z,shellNum)->OriginatingShellIds().size();
00510
00511 size_t index = 0;
00512
00513
00514
00515 while (provShellId != transitionManager->
00516 ReachableShell(Z,shellNum)->OriginatingShellId(index))
00517 {
00518 if(index == transitionSize-1)
00519 {
00520 break;
00521 }
00522 index++;
00523 }
00524
00525 G4double transitionEnergy = transitionManager->
00526 ReachableShell(Z,shellNum)->TransitionEnergy(index);
00527
00528 if (transitionEnergy < minGammaEnergy) return 0;
00529
00530
00531
00532 newShellId = transitionManager->
00533 ReachableShell(Z,shellNum)->OriginatingShellId(index);
00534
00535
00536 G4DynamicParticle* newPart = new G4DynamicParticle(G4Gamma::Gamma(),
00537 newGammaDirection,
00538 transitionEnergy);
00539 return newPart;
00540 }
00541
00542 G4DynamicParticle* G4UAtomicDeexcitation::GenerateAuger(G4int Z, G4int shellId)
00543 {
00544 if(!IsAugerActive()) {
00545
00546 return 0;
00547 }
00548
00549 if (shellId <=0 ) {
00550
00551
00552 G4Exception("G4UAtomicDeexcitation::GenerateAuger()","de0002",JustWarning, "Energy deposited locally");
00553
00554 return 0;
00555
00556 }
00557
00558
00559 G4int maxNumOfShells = transitionManager->NumberOfReachableAugerShells(Z);
00560
00561 const G4AugerTransition* refAugerTransition =
00562 transitionManager->ReachableAugerShell(Z,maxNumOfShells-1);
00563
00564
00565
00566
00567
00568 G4int shellNum = 0;
00569
00570 if ( shellId <= refAugerTransition->FinalShellId() )
00571
00572
00573 {
00574 G4int pippo = transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId();
00575 if (shellId != pippo ) {
00576 do {
00577 shellNum++;
00578 if(shellNum == maxNumOfShells)
00579 {
00580
00581 return 0;
00582 }
00583 }
00584 while (shellId != (transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId()) ) ;
00585 }
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595 G4int transitionLoopShellIndex = 0;
00596 G4double partSum = 0;
00597 const G4AugerTransition* anAugerTransition =
00598 transitionManager->ReachableAugerShell(Z,shellNum);
00599
00600
00601
00602
00603 G4int transitionSize =
00604 (anAugerTransition->TransitionOriginatingShellIds())->size();
00605 while (transitionLoopShellIndex < transitionSize) {
00606
00607 std::vector<G4int>::const_iterator pos =
00608 anAugerTransition->TransitionOriginatingShellIds()->begin();
00609
00610 G4int transitionLoopShellId = *(pos+transitionLoopShellIndex);
00611 G4int numberOfPossibleAuger =
00612 (anAugerTransition->AugerTransitionProbabilities(transitionLoopShellId))->size();
00613 G4int augerIndex = 0;
00614
00615
00616
00617 if (augerIndex < numberOfPossibleAuger) {
00618
00619 do
00620 {
00621 G4double thisProb = anAugerTransition->AugerTransitionProbability(augerIndex,
00622 transitionLoopShellId);
00623 partSum += thisProb;
00624 augerIndex++;
00625
00626 } while (augerIndex < numberOfPossibleAuger);
00627 }
00628 transitionLoopShellIndex++;
00629 }
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00655
00656 G4double totalVacancyAugerProbability = partSum;
00657
00658
00659
00660 G4int transitionRandomShellIndex = 0;
00661 G4int transitionRandomShellId = 1;
00662 G4int augerIndex = 0;
00663 partSum = 0;
00664 G4double partialProb = G4UniformRand();
00665
00666
00667 G4int numberOfPossibleAuger = 0;
00668
00669 G4bool foundFlag = false;
00670
00671 while (transitionRandomShellIndex < transitionSize) {
00672
00673 std::vector<G4int>::const_iterator pos =
00674 anAugerTransition->TransitionOriginatingShellIds()->begin();
00675
00676 transitionRandomShellId = *(pos+transitionRandomShellIndex);
00677
00678 augerIndex = 0;
00679 numberOfPossibleAuger = (anAugerTransition->
00680 AugerTransitionProbabilities(transitionRandomShellId))->size();
00681
00682 while (augerIndex < numberOfPossibleAuger) {
00683 G4double thisProb =anAugerTransition->AugerTransitionProbability(augerIndex,
00684 transitionRandomShellId);
00685
00686 partSum += thisProb;
00687
00688 if (partSum >= (partialProb*totalVacancyAugerProbability) ) {
00689 foundFlag = true;
00690 break;
00691 }
00692 augerIndex++;
00693 }
00694 if (partSum >= (partialProb*totalVacancyAugerProbability) ) {break;}
00695 transitionRandomShellIndex++;
00696 }
00697
00698
00699
00700
00701
00702 if (!foundFlag) {
00703
00704 return 0;
00705
00706 }
00707
00708
00709 G4double newcosTh = 1.-2.*G4UniformRand();
00710 G4double newsinTh = std::sqrt(1.-newcosTh*newcosTh);
00711 G4double newPhi = twopi*G4UniformRand();
00712
00713 G4double xDir = newsinTh*std::sin(newPhi);
00714 G4double yDir = newsinTh*std::cos(newPhi);
00715 G4double zDir = newcosTh;
00716
00717 G4ThreeVector newElectronDirection(xDir,yDir,zDir);
00718
00719
00720
00721
00722 G4double transitionEnergy = anAugerTransition->AugerTransitionEnergy(augerIndex, transitionRandomShellId);
00723
00724
00725
00726
00727
00728
00729 if (transitionEnergy < minElectronEnergy) {
00730
00731
00732
00733 return 0;
00734 }
00735
00736
00737
00738 newShellId = transitionRandomShellId;
00739
00740 return new G4DynamicParticle(G4Electron::Electron(),
00741 newElectronDirection,
00742 transitionEnergy);
00743 }
00744 else
00745 {
00746
00747
00748 return 0;
00749 }
00750 }