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 // $Id: G4PEEffectFluoModel.cc 73607 2013-09-02 10:04:03Z gcosmo $ 26 // 27 // 27 // ------------------------------------------- 28 // ------------------------------------------------------------------- 28 // 29 // 29 // GEANT4 Class file 30 // GEANT4 Class file 30 // 31 // 31 // 32 // 32 // File name: G4PEEffectFluoModel 33 // File name: G4PEEffectFluoModel 33 // 34 // 34 // Author: Vladimir Ivanchenko on base 35 // Author: Vladimir Ivanchenko on base of G4PEEffectModel 35 // 36 // 36 // Creation date: 13.06.2010 37 // Creation date: 13.06.2010 37 // 38 // 38 // Modifications: 39 // Modifications: 39 // 40 // 40 // Class Description: 41 // Class Description: 41 // Implementation of the photo-electric effect 42 // Implementation of the photo-electric effect with deexcitation 42 // 43 // 43 // ------------------------------------------- 44 // ------------------------------------------------------------------- 44 // 45 // 45 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 46 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 47 48 48 #include "G4PEEffectFluoModel.hh" 49 #include "G4PEEffectFluoModel.hh" 49 #include "G4PhysicalConstants.hh" 50 #include "G4PhysicalConstants.hh" 50 #include "G4SystemOfUnits.hh" 51 #include "G4SystemOfUnits.hh" 51 #include "G4Electron.hh" 52 #include "G4Electron.hh" 52 #include "G4Gamma.hh" 53 #include "G4Gamma.hh" 53 #include "Randomize.hh" 54 #include "Randomize.hh" 54 #include "G4Material.hh" 55 #include "G4Material.hh" 55 #include "G4DataVector.hh" 56 #include "G4DataVector.hh" 56 #include "G4ParticleChangeForGamma.hh" 57 #include "G4ParticleChangeForGamma.hh" 57 #include "G4VAtomDeexcitation.hh" 58 #include "G4VAtomDeexcitation.hh" 58 #include "G4LossTableManager.hh" 59 #include "G4LossTableManager.hh" 59 #include "G4SauterGavrilaAngularDistribution.h 60 #include "G4SauterGavrilaAngularDistribution.hh" 60 61 61 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 63 63 using namespace std; 64 using namespace std; 64 65 65 G4PEEffectFluoModel::G4PEEffectFluoModel(const 66 G4PEEffectFluoModel::G4PEEffectFluoModel(const G4String& nam) 66 : G4VEmModel(nam) 67 : G4VEmModel(nam) 67 { 68 { 68 theGamma = G4Gamma::Gamma(); 69 theGamma = G4Gamma::Gamma(); 69 theElectron = G4Electron::Electron(); 70 theElectron = G4Electron::Electron(); 70 fminimalEnergy = 1.0*CLHEP::eV; << 71 fminimalEnergy = 1.0*eV; 71 SetDeexcitationFlag(true); 72 SetDeexcitationFlag(true); >> 73 fParticleChange = 0; >> 74 fAtomDeexcitation = 0; 72 75 73 fSandiaCof.resize(4,0.0); 76 fSandiaCof.resize(4,0.0); 74 77 75 // default generator 78 // default generator 76 SetAngularDistribution(new G4SauterGavrilaAn 79 SetAngularDistribution(new G4SauterGavrilaAngularDistribution()); 77 } 80 } 78 81 79 //....oooOO0OOooo........oooOO0OOooo........oo 82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 80 83 81 G4PEEffectFluoModel::~G4PEEffectFluoModel() = << 84 G4PEEffectFluoModel::~G4PEEffectFluoModel() >> 85 {} 82 86 83 //....oooOO0OOooo........oooOO0OOooo........oo 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 84 88 85 void G4PEEffectFluoModel::Initialise(const G4P 89 void G4PEEffectFluoModel::Initialise(const G4ParticleDefinition*, 86 const G4DataVector&) 90 const G4DataVector&) 87 { 91 { 88 fAtomDeexcitation = G4LossTableManager::Inst 92 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation(); 89 fPEBelowKShell = G4EmParameters::Instance()- << 93 if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); } 90 if(nullptr == fParticleChange) { << 91 fParticleChange = GetParticleChangeForGamm << 92 } << 93 std::size_t nmat = G4Material::GetNumberOfMa << 94 fMatEnergyTh.resize(nmat, 0.0); << 95 for(std::size_t i=0; i<nmat; ++i) { << 96 fMatEnergyTh[i] = (*(G4Material::GetMateri << 97 ->GetSandiaTable()->GetSandiaCofForMater << 98 //G4cout << "G4PEEffectFluoModel::Initiali << 99 // << fMatEnergyTh[i]/eV << G4endl; << 100 } << 101 } 94 } 102 95 103 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo..... 104 97 105 G4double 98 G4double 106 G4PEEffectFluoModel::ComputeCrossSectionPerAto 99 G4PEEffectFluoModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, 107 G4double energy, 100 G4double energy, 108 G4double Z, G4double, 101 G4double Z, G4double, 109 G4double, G4double) 102 G4double, G4double) 110 { 103 { 111 // This method may be used only if G4Materia 104 // This method may be used only if G4MaterialCutsCouple pointer 112 // has been set properly 105 // has been set properly >> 106 113 CurrentCouple()->GetMaterial() 107 CurrentCouple()->GetMaterial() 114 ->GetSandiaTable()->GetSandiaCofPerAtom((G 108 ->GetSandiaTable()->GetSandiaCofPerAtom((G4int)Z, energy, fSandiaCof); 115 109 116 G4double x1 = 1 / energy; << 110 G4double energy2 = energy*energy; >> 111 G4double energy3 = energy*energy2; >> 112 G4double energy4 = energy2*energy2; 117 113 118 return x1 * (fSandiaCof[0] + x1 * (fSandiaCo << 114 return fSandiaCof[0]/energy + fSandiaCof[1]/energy2 + 119 x1 * (fSandiaCof[2] + x1 * fSandiaCof[3])) << 115 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4; 120 } 116 } 121 117 122 //....oooOO0OOooo........oooOO0OOooo........oo 118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 119 124 G4double 120 G4double 125 G4PEEffectFluoModel::CrossSectionPerVolume(con 121 G4PEEffectFluoModel::CrossSectionPerVolume(const G4Material* material, 126 const G4ParticleDefinition*, 122 const G4ParticleDefinition*, 127 G4double energy, 123 G4double energy, 128 G4double, G4double) 124 G4double, G4double) 129 { 125 { 130 // This method may be used only if G4Materia << 131 // has been set properly << 132 energy = std::max(energy, fMatEnergyTh[mater << 133 const G4double* SandiaCof = 126 const G4double* SandiaCof = 134 material->GetSandiaTable()->GetSandiaCofFo 127 material->GetSandiaTable()->GetSandiaCofForMaterial(energy); 135 << 128 136 G4double x1 = 1 / energy; << 129 G4double energy2 = energy*energy; 137 << 130 G4double energy3 = energy*energy2; 138 return x1 * (SandiaCof[0] + x1 * (SandiaCof[ << 131 G4double energy4 = energy2*energy2; 139 x1 * (SandiaCof[2] + x1 * SandiaCof[3]))); << 132 >> 133 return SandiaCof[0]/energy + SandiaCof[1]/energy2 + >> 134 SandiaCof[2]/energy3 + SandiaCof[3]/energy4; 140 } 135 } 141 136 142 //....oooOO0OOooo........oooOO0OOooo........oo 137 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 143 138 144 void 139 void 145 G4PEEffectFluoModel::SampleSecondaries(std::ve 140 G4PEEffectFluoModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 146 const G4MaterialCutsCouple* cou 141 const G4MaterialCutsCouple* couple, 147 const G4DynamicParticle* aDynam 142 const G4DynamicParticle* aDynamicPhoton, 148 G4double, 143 G4double, 149 G4double) 144 G4double) 150 { 145 { 151 SetCurrentCouple(couple); 146 SetCurrentCouple(couple); 152 const G4Material* aMaterial = couple->GetMat 147 const G4Material* aMaterial = couple->GetMaterial(); 153 148 154 G4double energy = aDynamicPhoton->GetKinetic 149 G4double energy = aDynamicPhoton->GetKineticEnergy(); 155 150 156 // select randomly one element constituing t 151 // select randomly one element constituing the material. 157 const G4Element* anElement = SelectRandomAto 152 const G4Element* anElement = SelectRandomAtom(aMaterial,theGamma,energy); 158 153 159 // 154 // 160 // Photo electron 155 // Photo electron 161 // 156 // 162 157 163 // Select atomic shell 158 // Select atomic shell 164 G4int nShells = anElement->GetNbOfAtomicShel 159 G4int nShells = anElement->GetNbOfAtomicShells(); 165 G4int i = 0; 160 G4int i = 0; 166 for(; i<nShells; ++i) { 161 for(; i<nShells; ++i) { 167 /* 162 /* 168 G4cout << "i= " << i << " E(eV)= " << ener 163 G4cout << "i= " << i << " E(eV)= " << energy/eV 169 << " Eb(eV)= " << anElement->GetAtomi 164 << " Eb(eV)= " << anElement->GetAtomicShell(i)/eV 170 << " " << anElement->GetName() 165 << " " << anElement->GetName() 171 << G4endl; 166 << G4endl; 172 */ 167 */ 173 if(energy >= anElement->GetAtomicShell(i)) 168 if(energy >= anElement->GetAtomicShell(i)) { break; } 174 } 169 } 175 170 176 G4double edep = energy; 171 G4double edep = energy; 177 172 178 // photo-electron is not sampled if shell is << 173 // Normally one shell is available 179 // the flag of photoeffect is "false" and sh << 174 if (i < nShells) { 180 if ( (fPEBelowKShell || 0 == i) && i < nShel << 175 181 << 182 G4double bindingEnergy = anElement->GetAto 176 G4double bindingEnergy = anElement->GetAtomicShell(i); 183 edep = bindingEnergy; 177 edep = bindingEnergy; 184 G4double esec = 0.0; 178 G4double esec = 0.0; 185 179 186 // sample deexcitation cascade << 180 // sample deexcitation 187 // 181 // 188 if(nullptr != fAtomDeexcitation) { << 182 if(fAtomDeexcitation) { 189 G4int index = couple->GetIndex(); 183 G4int index = couple->GetIndex(); 190 if(fAtomDeexcitation->CheckDeexcitationA 184 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) { 191 G4int Z = G4lrint(anElement->GetZ()); 185 G4int Z = G4lrint(anElement->GetZ()); 192 auto as = (G4AtomicShellEnumerator)(i); << 186 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i); 193 const G4AtomicShell* shell = fAtomDeexcitati 187 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as); 194 G4double eshell = shell->BindingEnergy 188 G4double eshell = shell->BindingEnergy(); 195 if(eshell > bindingEnergy && eshell <= 189 if(eshell > bindingEnergy && eshell <= energy) { 196 bindingEnergy = eshell; 190 bindingEnergy = eshell; 197 edep = eshell; 191 edep = eshell; 198 } 192 } 199 std::size_t nbefore = fvect->size(); << 193 size_t nbefore = fvect->size(); 200 fAtomDeexcitation->GenerateParticles(fvect, 194 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index); 201 std::size_t nafter = fvect->size(); << 195 size_t nafter = fvect->size(); 202 for (std::size_t j=nbefore; j<nafter; ++j) { << 196 if(nafter > nbefore) { 203 G4double e = ((*fvect)[j])->GetKineticEner << 197 for (size_t j=nbefore; j<nafter; ++j) { 204 if(esec + e > edep) { << 198 G4double e = ((*fvect)[j])->GetKineticEnergy(); 205 // correct energy in order to have energ << 199 if(esec + e > edep) { 206 e = edep - esec; << 200 /* 207 ((*fvect)[j])->SetKineticEnergy(e); << 208 esec += e; << 209 /* << 210 G4cout << "### G4PEffectFluoModel Edep 201 G4cout << "### G4PEffectFluoModel Edep(eV)= " << edep/eV 211 << " Esec(eV)= " << esec/eV 202 << " Esec(eV)= " << esec/eV 212 << " E["<< j << "](eV)= " << e/eV 203 << " E["<< j << "](eV)= " << e/eV 213 << " N= " << nafter 204 << " N= " << nafter 214 << " Z= " << Z << " shell= " << i 205 << " Z= " << Z << " shell= " << i 215 << " Ebind(keV)= " << bindingEnergy/ 206 << " Ebind(keV)= " << bindingEnergy/keV 216 << " Eshell(keV)= " << shell->Bindin 207 << " Eshell(keV)= " << shell->BindingEnergy()/keV 217 << G4endl; 208 << G4endl; 218 */ << 209 */ 219 // delete the rest of secondaries (shoul << 210 for (size_t jj=j; jj<nafter; ++jj) { delete (*fvect)[jj]; } 220 for (std::size_t jj=nafter-1; jj>j; --jj << 211 for (size_t jj=j; jj<nafter; ++jj) { fvect->pop_back(); } 221 delete (*fvect)[jj]; << 212 break; 222 fvect->pop_back(); << 223 } 213 } 224 break; << 214 esec += e; 225 } << 215 } 226 esec += e; << 227 } 216 } 228 edep -= esec; 217 edep -= esec; 229 } 218 } 230 } 219 } 231 // create photo electron 220 // create photo electron 232 // 221 // 233 G4double elecKineEnergy = energy - binding 222 G4double elecKineEnergy = energy - bindingEnergy; 234 if (elecKineEnergy > fminimalEnergy) { 223 if (elecKineEnergy > fminimalEnergy) { 235 auto aParticle = new G4DynamicParticle(t << 224 G4DynamicParticle* aParticle = new G4DynamicParticle(theElectron, 236 GetAngularDistribution()->SampleDirection(aD 225 GetAngularDistribution()->SampleDirection(aDynamicPhoton, 237 elecKineEnergy, 226 elecKineEnergy, 238 i, couple->GetMaterial()), 227 i, couple->GetMaterial()), 239 elecKineEnergy); 228 elecKineEnergy); 240 fvect->push_back(aParticle); 229 fvect->push_back(aParticle); 241 } else { 230 } else { 242 edep += elecKineEnergy; 231 edep += elecKineEnergy; 243 elecKineEnergy = 0.0; 232 elecKineEnergy = 0.0; 244 } 233 } 245 if(std::abs(energy - elecKineEnergy - esec << 234 if(fabs(energy - elecKineEnergy - esec - edep) > eV) { 246 G4cout << "### G4PEffectFluoModel dE(eV) 235 G4cout << "### G4PEffectFluoModel dE(eV)= " 247 << (energy - elecKineEnergy - esec - ed 236 << (energy - elecKineEnergy - esec - edep)/eV 248 << " shell= " << i 237 << " shell= " << i 249 << " E(keV)= " << energy/keV 238 << " E(keV)= " << energy/keV 250 << " Ebind(keV)= " << bindingEnergy/ke 239 << " Ebind(keV)= " << bindingEnergy/keV 251 << " Ee(keV)= " << elecKineEnergy/keV 240 << " Ee(keV)= " << elecKineEnergy/keV 252 << " Esec(keV)= " << esec/keV 241 << " Esec(keV)= " << esec/keV 253 << " Edep(keV)= " << edep/keV 242 << " Edep(keV)= " << edep/keV 254 << G4endl; 243 << G4endl; 255 } 244 } 256 } 245 } 257 246 258 // kill primary photon 247 // kill primary photon 259 fParticleChange->SetProposedKineticEnergy(0. 248 fParticleChange->SetProposedKineticEnergy(0.); 260 fParticleChange->ProposeTrackStatus(fStopAnd 249 fParticleChange->ProposeTrackStatus(fStopAndKill); 261 if(edep > 0.0) { 250 if(edep > 0.0) { 262 fParticleChange->ProposeLocalEnergyDeposit 251 fParticleChange->ProposeLocalEnergyDeposit(edep); 263 } 252 } 264 } 253 } 265 254 266 //....oooOO0OOooo........oooOO0OOooo........oo 255 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 267 256