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 #include "G4LivermorePolarizedPhotoElectricModel.hh"
00033 #include "G4PhysicalConstants.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4LossTableManager.hh"
00036 #include "G4VAtomDeexcitation.hh"
00037 #include "G4PhotoElectricAngularGeneratorPolarized.hh"
00038 #include "G4SauterGavrilaAngularDistribution.hh"
00039 #include "G4AtomicShell.hh"
00040 #include "G4ProductionCutsTable.hh"
00041 #include "G4Gamma.hh"
00042
00043
00044
00045 using namespace std;
00046
00047
00048
00049 G4LivermorePolarizedPhotoElectricModel::G4LivermorePolarizedPhotoElectricModel(
00050 const G4String& nam)
00051 :G4VEmModel(nam),fParticleChange(0),
00052 crossSectionHandler(0), shellCrossSectionHandler(0)
00053 {
00054 lowEnergyLimit = 10 * eV;
00055 highEnergyLimit = 100 * GeV;
00056
00057 verboseLevel= 0;
00058
00059
00060
00061
00062
00063
00064
00065 theGamma = G4Gamma::Gamma();
00066 theElectron = G4Electron::Electron();
00067
00068
00069 SetAngularDistribution(new G4SauterGavrilaAngularDistribution());
00070
00071
00072
00073 SetDeexcitationFlag(true);
00074 fAtomDeexcitation = 0;
00075 fDeexcitationActive = false;
00076
00077 if (verboseLevel > 1) {
00078 G4cout << "Livermore Polarized PhotoElectric is constructed " << G4endl
00079 << "Energy range: "
00080 << lowEnergyLimit / keV << " keV - "
00081 << highEnergyLimit / GeV << " GeV"
00082 << G4endl;
00083 }
00084 }
00085
00086
00087
00088 G4LivermorePolarizedPhotoElectricModel::~G4LivermorePolarizedPhotoElectricModel()
00089 {
00090 delete crossSectionHandler;
00091 delete shellCrossSectionHandler;
00092 }
00093
00094
00095
00096 void
00097 G4LivermorePolarizedPhotoElectricModel::Initialise(
00098 const G4ParticleDefinition*,
00099 const G4DataVector&)
00100 {
00101 if (verboseLevel > 3) {
00102 G4cout << "Calling G4LivermorePolarizedPhotoElectricModel::Initialise()" << G4endl;
00103 }
00104 if (crossSectionHandler)
00105 {
00106 crossSectionHandler->Clear();
00107 delete crossSectionHandler;
00108 }
00109
00110 if (shellCrossSectionHandler)
00111 {
00112 shellCrossSectionHandler->Clear();
00113 delete shellCrossSectionHandler;
00114 }
00115
00116
00117
00118 crossSectionHandler = new G4CrossSectionHandler;
00119 crossSectionHandler->Clear();
00120 G4String crossSectionFile = "phot/pe-cs-";
00121 crossSectionHandler->LoadData(crossSectionFile);
00122
00123 shellCrossSectionHandler = new G4CrossSectionHandler();
00124 shellCrossSectionHandler->Clear();
00125 G4String shellCrossSectionFile = "phot/pe-ss-cs-";
00126 shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
00127
00128 fParticleChange = GetParticleChangeForGamma();
00129
00130 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00131 if(fAtomDeexcitation) {
00132 fDeexcitationActive = fAtomDeexcitation->IsFluoActive();
00133 }
00134 if (verboseLevel > 1) {
00135 G4cout << "Livermore Polarized PhotoElectric model is initialized "
00136 << G4endl
00137 << "Energy range: "
00138 << LowEnergyLimit() / keV << " keV - "
00139 << HighEnergyLimit() / GeV << " GeV"
00140 << G4endl;
00141 }
00142 }
00143
00144
00145
00146 G4double G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom(
00147 const G4ParticleDefinition*,
00148 G4double GammaEnergy,
00149 G4double Z, G4double,
00150 G4double, G4double)
00151 {
00152 G4double cs = crossSectionHandler->FindValue(G4lrint(Z), GammaEnergy);
00153
00154 if (verboseLevel > 3) {
00155 G4cout << "G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom()"
00156 << G4endl;
00157 G4cout << " E(keV)= " << GammaEnergy/keV << " Z= " << Z
00158 << " CrossSection(barn)= "
00159 << cs/barn << G4endl;
00160 }
00161 return cs;
00162 }
00163
00164
00165
00166 void G4LivermorePolarizedPhotoElectricModel::SampleSecondaries(
00167 std::vector<G4DynamicParticle*>* fvect,
00168 const G4MaterialCutsCouple* couple,
00169 const G4DynamicParticle* aDynamicGamma,
00170 G4double,
00171 G4double)
00172 {
00173 if (verboseLevel > 3) {
00174 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricModel"
00175 << G4endl;
00176 }
00177
00178 G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
00179
00180
00181 fParticleChange->SetProposedKineticEnergy(0.);
00182 fParticleChange->ProposeTrackStatus(fStopAndKill);
00183
00184
00185
00186 if (photonEnergy <= lowEnergyLimit)
00187 {
00188 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00189 return;
00190 }
00191
00192
00193
00194 const G4Element* elm =
00195 SelectRandomAtom(couple->GetMaterial(),theGamma,photonEnergy);
00196 G4int Z = G4lrint(elm->GetZ());
00197
00198
00199
00200 size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
00201
00202
00203 G4double bindingEnergy;
00204 const G4AtomicShell* shell = 0;
00205 if(fDeexcitationActive) {
00206 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIndex);
00207 shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00208 bindingEnergy = shell->BindingEnergy();
00209 } else {
00210 G4int nshells = elm->GetNbOfAtomicShells() - 1;
00211 if(G4int(shellIndex) > nshells) { shellIndex = std::max(0, nshells); }
00212 bindingEnergy = elm->GetAtomicShell(shellIndex);
00213 }
00214
00215
00216
00217
00218 if(photonEnergy < bindingEnergy) {
00219 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
00220 return;
00221 }
00222
00223
00224
00225
00226
00227 G4double eKineticEnergy = photonEnergy - bindingEnergy;
00228 G4double edep = bindingEnergy;
00229
00230
00231 G4ThreeVector electronDirection =
00232 GetAngularDistribution()->SampleDirection(aDynamicGamma,
00233 eKineticEnergy,
00234 shellIndex,
00235 couple->GetMaterial());
00236
00237
00238 G4DynamicParticle* electron = new G4DynamicParticle (theElectron,
00239 electronDirection,
00240 eKineticEnergy);
00241 fvect->push_back(electron);
00242
00243
00244 if(fDeexcitationActive) {
00245 G4int index = couple->GetIndex();
00246 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00247 size_t nbefore = fvect->size();
00248
00249 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00250 size_t nafter = fvect->size();
00251 if(nafter > nbefore) {
00252 for (size_t j=nbefore; j<nafter; ++j) {
00253 edep -= ((*fvect)[j])->GetKineticEnergy();
00254 }
00255 }
00256 }
00257 }
00258
00259
00260 if(edep > 0.0) {
00261 fParticleChange->ProposeLocalEnergyDeposit(edep);
00262 }
00263 }
00264
00265