Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4PEEffectFluoModel 32 // File name: G4PEEffectFluoModel 33 // 33 // 34 // Author: Vladimir Ivanchenko on base 34 // Author: Vladimir Ivanchenko on base of G4PEEffectModel 35 // 35 // 36 // Creation date: 13.06.2010 36 // Creation date: 13.06.2010 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // 40 // Class Description: 40 // Class Description: 41 // Implementation of the photo-electric effect 41 // Implementation of the photo-electric effect with deexcitation 42 // 42 // 43 // ------------------------------------------- 43 // ------------------------------------------------------------------- 44 // 44 // 45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 47 48 #include "G4PEEffectFluoModel.hh" 48 #include "G4PEEffectFluoModel.hh" 49 #include "G4PhysicalConstants.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "G4Electron.hh" 51 #include "G4Electron.hh" 52 #include "G4Gamma.hh" 52 #include "G4Gamma.hh" 53 #include "Randomize.hh" 53 #include "Randomize.hh" 54 #include "G4Material.hh" 54 #include "G4Material.hh" 55 #include "G4DataVector.hh" 55 #include "G4DataVector.hh" 56 #include "G4ParticleChangeForGamma.hh" 56 #include "G4ParticleChangeForGamma.hh" 57 #include "G4VAtomDeexcitation.hh" 57 #include "G4VAtomDeexcitation.hh" 58 #include "G4LossTableManager.hh" 58 #include "G4LossTableManager.hh" 59 #include "G4SauterGavrilaAngularDistribution.h 59 #include "G4SauterGavrilaAngularDistribution.hh" 60 60 61 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 62 63 using namespace std; 63 using namespace std; 64 64 65 G4PEEffectFluoModel::G4PEEffectFluoModel(const 65 G4PEEffectFluoModel::G4PEEffectFluoModel(const G4String& nam) 66 : G4VEmModel(nam) 66 : G4VEmModel(nam) 67 { 67 { 68 theGamma = G4Gamma::Gamma(); 68 theGamma = G4Gamma::Gamma(); 69 theElectron = G4Electron::Electron(); 69 theElectron = G4Electron::Electron(); 70 fminimalEnergy = 1.0*CLHEP::eV; << 70 fminimalEnergy = 1.0*eV; 71 SetDeexcitationFlag(true); 71 SetDeexcitationFlag(true); >> 72 fParticleChange = nullptr; >> 73 fAtomDeexcitation = nullptr; 72 74 73 fSandiaCof.resize(4,0.0); 75 fSandiaCof.resize(4,0.0); 74 76 75 // default generator 77 // default generator 76 SetAngularDistribution(new G4SauterGavrilaAn 78 SetAngularDistribution(new G4SauterGavrilaAngularDistribution()); 77 } 79 } 78 80 79 //....oooOO0OOooo........oooOO0OOooo........oo 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 82 81 G4PEEffectFluoModel::~G4PEEffectFluoModel() = << 83 G4PEEffectFluoModel::~G4PEEffectFluoModel() >> 84 {} 82 85 83 //....oooOO0OOooo........oooOO0OOooo........oo 86 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 84 87 85 void G4PEEffectFluoModel::Initialise(const G4P 88 void G4PEEffectFluoModel::Initialise(const G4ParticleDefinition*, 86 const G4DataVector&) 89 const G4DataVector&) 87 { 90 { 88 fAtomDeexcitation = G4LossTableManager::Inst 91 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 89 fPEBelowKShell = G4EmParameters::Instance()- << 92 if(nullptr == fParticleChange) { fParticleChange = GetParticleChangeForGamma(); } 90 if(nullptr == fParticleChange) { << 93 size_t nmat = G4Material::GetNumberOfMaterials(); 91 fParticleChange = GetParticleChangeForGamm << 92 } << 93 std::size_t nmat = G4Material::GetNumberOfMa << 94 fMatEnergyTh.resize(nmat, 0.0); 94 fMatEnergyTh.resize(nmat, 0.0); 95 for(std::size_t i=0; i<nmat; ++i) { << 95 for(size_t i=0; i<nmat; ++i) { 96 fMatEnergyTh[i] = (*(G4Material::GetMateri 96 fMatEnergyTh[i] = (*(G4Material::GetMaterialTable()))[i] 97 ->GetSandiaTable()->GetSandiaCofForMater 97 ->GetSandiaTable()->GetSandiaCofForMaterial(0, 0); 98 //G4cout << "G4PEEffectFluoModel::Initiali 98 //G4cout << "G4PEEffectFluoModel::Initialise Eth(eV)= " 99 // << fMatEnergyTh[i]/eV << G4endl; 99 // << fMatEnergyTh[i]/eV << G4endl; 100 } 100 } 101 } 101 } 102 102 103 //....oooOO0OOooo........oooOO0OOooo........oo 103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 104 104 105 G4double 105 G4double 106 G4PEEffectFluoModel::ComputeCrossSectionPerAto 106 G4PEEffectFluoModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 107 G4double energy, 107 G4double energy, 108 G4double Z, G4double, 108 G4double Z, G4double, 109 G4double, G4double) 109 G4double, G4double) 110 { 110 { 111 // This method may be used only if G4Materia 111 // This method may be used only if G4MaterialCutsCouple pointer 112 // has been set properly 112 // has been set properly 113 CurrentCouple()->GetMaterial() 113 CurrentCouple()->GetMaterial() 114 ->GetSandiaTable()->GetSandiaCofPerAtom((G 114 ->GetSandiaTable()->GetSandiaCofPerAtom((G4int)Z, energy, fSandiaCof); 115 115 116 G4double x1 = 1 / energy; << 116 G4double energy2 = energy*energy; >> 117 G4double energy3 = energy*energy2; >> 118 G4double energy4 = energy2*energy2; 117 119 118 return x1 * (fSandiaCof[0] + x1 * (fSandiaCo << 120 return fSandiaCof[0]/energy + fSandiaCof[1]/energy2 + 119 x1 * (fSandiaCof[2] + x1 * fSandiaCof[3])) << 121 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4; 120 } 122 } 121 123 122 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 125 124 G4double 126 G4double 125 G4PEEffectFluoModel::CrossSectionPerVolume(con 127 G4PEEffectFluoModel::CrossSectionPerVolume(const G4Material* material, 126 const G4ParticleDefinition*, 128 const G4ParticleDefinition*, 127 G4double energy, 129 G4double energy, 128 G4double, G4double) 130 G4double, G4double) 129 { 131 { 130 // This method may be used only if G4Materia 132 // This method may be used only if G4MaterialCutsCouple pointer 131 // has been set properly 133 // has been set properly 132 energy = std::max(energy, fMatEnergyTh[mater 134 energy = std::max(energy, fMatEnergyTh[material->GetIndex()]); 133 const G4double* SandiaCof = 135 const G4double* SandiaCof = 134 material->GetSandiaTable()->GetSandiaCofFo 136 material->GetSandiaTable()->GetSandiaCofForMaterial(energy); 135 << 137 136 G4double x1 = 1 / energy; << 138 G4double energy2 = energy*energy; 137 << 139 G4double energy3 = energy*energy2; 138 return x1 * (SandiaCof[0] + x1 * (SandiaCof[ << 140 G4double energy4 = energy2*energy2; 139 x1 * (SandiaCof[2] + x1 * SandiaCof[3]))); << 141 >> 142 return SandiaCof[0]/energy + SandiaCof[1]/energy2 + >> 143 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 140 } 144 } 141 145 142 //....oooOO0OOooo........oooOO0OOooo........oo 146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 143 147 144 void 148 void 145 G4PEEffectFluoModel::SampleSecondaries(std::ve 149 G4PEEffectFluoModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 146 const G4MaterialCutsCouple* cou 150 const G4MaterialCutsCouple* couple, 147 const G4DynamicParticle* aDynam 151 const G4DynamicParticle* aDynamicPhoton, 148 G4double, 152 G4double, 149 G4double) 153 G4double) 150 { 154 { 151 SetCurrentCouple(couple); 155 SetCurrentCouple(couple); 152 const G4Material* aMaterial = couple->GetMat 156 const G4Material* aMaterial = couple->GetMaterial(); 153 157 154 G4double energy = aDynamicPhoton->GetKinetic 158 G4double energy = aDynamicPhoton->GetKineticEnergy(); 155 159 156 // select randomly one element constituing t 160 // select randomly one element constituing the material. 157 const G4Element* anElement = SelectRandomAto 161 const G4Element* anElement = SelectRandomAtom(aMaterial,theGamma,energy); 158 162 159 // 163 // 160 // Photo electron 164 // Photo electron 161 // 165 // 162 166 163 // Select atomic shell 167 // Select atomic shell 164 G4int nShells = anElement->GetNbOfAtomicShel 168 G4int nShells = anElement->GetNbOfAtomicShells(); 165 G4int i = 0; 169 G4int i = 0; 166 for(; i<nShells; ++i) { 170 for(; i<nShells; ++i) { 167 /* 171 /* 168 G4cout << "i= " << i << " E(eV)= " << ener 172 G4cout << "i= " << i << " E(eV)= " << energy/eV 169 << " Eb(eV)= " << anElement->GetAtomi 173 << " Eb(eV)= " << anElement->GetAtomicShell(i)/eV 170 << " " << anElement->GetName() 174 << " " << anElement->GetName() 171 << G4endl; 175 << G4endl; 172 */ 176 */ 173 if(energy >= anElement->GetAtomicShell(i)) 177 if(energy >= anElement->GetAtomicShell(i)) { break; } 174 } 178 } 175 179 176 G4double edep = energy; 180 G4double edep = energy; 177 181 178 // photo-electron is not sampled if shell is << 182 // Normally one shell is available 179 // the flag of photoeffect is "false" and sh << 183 if (i < nShells) { 180 if ( (fPEBelowKShell || 0 == i) && i < nShel << 184 181 << 182 G4double bindingEnergy = anElement->GetAto 185 G4double bindingEnergy = anElement->GetAtomicShell(i); 183 edep = bindingEnergy; 186 edep = bindingEnergy; 184 G4double esec = 0.0; 187 G4double esec = 0.0; 185 188 186 // sample deexcitation cascade << 189 // sample deexcitation 187 // 190 // 188 if(nullptr != fAtomDeexcitation) { << 191 if(fAtomDeexcitation) { 189 G4int index = couple->GetIndex(); 192 G4int index = couple->GetIndex(); 190 if(fAtomDeexcitation->CheckDeexcitationA 193 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 191 G4int Z = G4lrint(anElement->GetZ()); 194 G4int Z = G4lrint(anElement->GetZ()); 192 auto as = (G4AtomicShellEnumerator)(i); << 195 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i); 193 const G4AtomicShell* shell = fAtomDeexcitati 196 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 194 G4double eshell = shell->BindingEnergy 197 G4double eshell = shell->BindingEnergy(); 195 if(eshell > bindingEnergy && eshell <= 198 if(eshell > bindingEnergy && eshell <= energy) { 196 bindingEnergy = eshell; 199 bindingEnergy = eshell; 197 edep = eshell; 200 edep = eshell; 198 } 201 } 199 std::size_t nbefore = fvect->size(); << 202 G4int nbefore = fvect->size(); 200 fAtomDeexcitation->GenerateParticles(fvect, 203 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 201 std::size_t nafter = fvect->size(); << 204 G4int nafter = fvect->size(); 202 for (std::size_t j=nbefore; j<nafter; ++j) { << 205 for (G4int j=nbefore; j<nafter; ++j) { 203 G4double e = ((*fvect)[j])->GetKineticEner 206 G4double e = ((*fvect)[j])->GetKineticEnergy(); 204 if(esec + e > edep) { 207 if(esec + e > edep) { 205 // correct energy in order to have energ 208 // correct energy in order to have energy balance 206 e = edep - esec; 209 e = edep - esec; 207 ((*fvect)[j])->SetKineticEnergy(e); 210 ((*fvect)[j])->SetKineticEnergy(e); 208 esec += e; 211 esec += e; 209 /* 212 /* 210 G4cout << "### G4PEffectFluoModel Edep 213 G4cout << "### G4PEffectFluoModel Edep(eV)= " << edep/eV 211 << " Esec(eV)= " << esec/eV 214 << " Esec(eV)= " << esec/eV 212 << " E["<< j << "](eV)= " << e/eV 215 << " E["<< j << "](eV)= " << e/eV 213 << " N= " << nafter 216 << " N= " << nafter 214 << " Z= " << Z << " shell= " << i 217 << " Z= " << Z << " shell= " << i 215 << " Ebind(keV)= " << bindingEnergy/ 218 << " Ebind(keV)= " << bindingEnergy/keV 216 << " Eshell(keV)= " << shell->Bindin 219 << " Eshell(keV)= " << shell->BindingEnergy()/keV 217 << G4endl; 220 << G4endl; 218 */ 221 */ 219 // delete the rest of secondaries (shoul 222 // delete the rest of secondaries (should not happens) 220 for (std::size_t jj=nafter-1; jj>j; --jj << 223 for (G4int jj=nafter-1; jj>j; --jj) { 221 delete (*fvect)[jj]; 224 delete (*fvect)[jj]; 222 fvect->pop_back(); 225 fvect->pop_back(); 223 } 226 } 224 break; 227 break; 225 } 228 } 226 esec += e; 229 esec += e; 227 } 230 } 228 edep -= esec; 231 edep -= esec; 229 } 232 } 230 } 233 } 231 // create photo electron 234 // create photo electron 232 // 235 // 233 G4double elecKineEnergy = energy - binding 236 G4double elecKineEnergy = energy - bindingEnergy; 234 if (elecKineEnergy > fminimalEnergy) { 237 if (elecKineEnergy > fminimalEnergy) { 235 auto aParticle = new G4DynamicParticle(t << 238 G4DynamicParticle* aParticle = new G4DynamicParticle(theElectron, 236 GetAngularDistribution()->SampleDirection(aD 239 GetAngularDistribution()->SampleDirection(aDynamicPhoton, 237 elecKineEnergy, 240 elecKineEnergy, 238 i, couple->GetMaterial()), 241 i, couple->GetMaterial()), 239 elecKineEnergy); 242 elecKineEnergy); 240 fvect->push_back(aParticle); 243 fvect->push_back(aParticle); 241 } else { 244 } else { 242 edep += elecKineEnergy; 245 edep += elecKineEnergy; 243 elecKineEnergy = 0.0; 246 elecKineEnergy = 0.0; 244 } 247 } 245 if(std::abs(energy - elecKineEnergy - esec 248 if(std::abs(energy - elecKineEnergy - esec - edep) > CLHEP::eV) { 246 G4cout << "### G4PEffectFluoModel dE(eV) 249 G4cout << "### G4PEffectFluoModel dE(eV)= " 247 << (energy - elecKineEnergy - esec - ed 250 << (energy - elecKineEnergy - esec - edep)/eV 248 << " shell= " << i 251 << " shell= " << i 249 << " E(keV)= " << energy/keV 252 << " E(keV)= " << energy/keV 250 << " Ebind(keV)= " << bindingEnergy/ke 253 << " Ebind(keV)= " << bindingEnergy/keV 251 << " Ee(keV)= " << elecKineEnergy/keV 254 << " Ee(keV)= " << elecKineEnergy/keV 252 << " Esec(keV)= " << esec/keV 255 << " Esec(keV)= " << esec/keV 253 << " Edep(keV)= " << edep/keV 256 << " Edep(keV)= " << edep/keV 254 << G4endl; 257 << G4endl; 255 } 258 } 256 } 259 } 257 260 258 // kill primary photon 261 // kill primary photon 259 fParticleChange->SetProposedKineticEnergy(0. 262 fParticleChange->SetProposedKineticEnergy(0.); 260 fParticleChange->ProposeTrackStatus(fStopAnd 263 fParticleChange->ProposeTrackStatus(fStopAndKill); 261 if(edep > 0.0) { 264 if(edep > 0.0) { 262 fParticleChange->ProposeLocalEnergyDeposit 265 fParticleChange->ProposeLocalEnergyDeposit(edep); 263 } 266 } 264 } 267 } 265 268 266 //....oooOO0OOooo........oooOO0OOooo........oo 269 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 267 270