| 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 106442 2017-10-10 08:00:30Z 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 = nullptr;
>> 74 fAtomDeexcitation = nullptr;
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(nullptr == fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
90 if(nullptr == fParticleChange) { << 94 size_t nmat = G4Material::GetNumberOfMaterials();
91 fParticleChange = GetParticleChangeForGamm <<
92 } <<
93 std::size_t nmat = G4Material::GetNumberOfMa <<
94 fMatEnergyTh.resize(nmat, 0.0); 95 fMatEnergyTh.resize(nmat, 0.0);
95 for(std::size_t i=0; i<nmat; ++i) { << 96 for(size_t i=0; i<nmat; ++i) {
96 fMatEnergyTh[i] = (*(G4Material::GetMateri 97 fMatEnergyTh[i] = (*(G4Material::GetMaterialTable()))[i]
97 ->GetSandiaTable()->GetSandiaCofForMater 98 ->GetSandiaTable()->GetSandiaCofForMaterial(0, 0);
98 //G4cout << "G4PEEffectFluoModel::Initiali 99 //G4cout << "G4PEEffectFluoModel::Initialise Eth(eV)= "
99 // << fMatEnergyTh[i]/eV << G4endl; 100 // << fMatEnergyTh[i]/eV << G4endl;
100 } 101 }
101 } 102 }
102 103
103 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
104 105
105 G4double 106 G4double
106 G4PEEffectFluoModel::ComputeCrossSectionPerAto 107 G4PEEffectFluoModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
107 G4double energy, 108 G4double energy,
108 G4double Z, G4double, 109 G4double Z, G4double,
109 G4double, G4double) 110 G4double, G4double)
110 { 111 {
111 // This method may be used only if G4Materia 112 // This method may be used only if G4MaterialCutsCouple pointer
112 // has been set properly 113 // has been set properly
113 CurrentCouple()->GetMaterial() 114 CurrentCouple()->GetMaterial()
114 ->GetSandiaTable()->GetSandiaCofPerAtom((G 115 ->GetSandiaTable()->GetSandiaCofPerAtom((G4int)Z, energy, fSandiaCof);
115 116
116 G4double x1 = 1 / energy; << 117 G4double energy2 = energy*energy;
>> 118 G4double energy3 = energy*energy2;
>> 119 G4double energy4 = energy2*energy2;
117 120
118 return x1 * (fSandiaCof[0] + x1 * (fSandiaCo << 121 return fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
119 x1 * (fSandiaCof[2] + x1 * fSandiaCof[3])) << 122 fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4;
120 } 123 }
121 124
122 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
123 126
124 G4double 127 G4double
125 G4PEEffectFluoModel::CrossSectionPerVolume(con 128 G4PEEffectFluoModel::CrossSectionPerVolume(const G4Material* material,
126 const G4ParticleDefinition*, 129 const G4ParticleDefinition*,
127 G4double energy, 130 G4double energy,
128 G4double, G4double) 131 G4double, G4double)
129 { 132 {
130 // This method may be used only if G4Materia 133 // This method may be used only if G4MaterialCutsCouple pointer
131 // has been set properly 134 // has been set properly
132 energy = std::max(energy, fMatEnergyTh[mater 135 energy = std::max(energy, fMatEnergyTh[material->GetIndex()]);
133 const G4double* SandiaCof = 136 const G4double* SandiaCof =
134 material->GetSandiaTable()->GetSandiaCofFo 137 material->GetSandiaTable()->GetSandiaCofForMaterial(energy);
135 << 138
136 G4double x1 = 1 / energy; << 139 G4double energy2 = energy*energy;
137 << 140 G4double energy3 = energy*energy2;
138 return x1 * (SandiaCof[0] + x1 * (SandiaCof[ << 141 G4double energy4 = energy2*energy2;
139 x1 * (SandiaCof[2] + x1 * SandiaCof[3]))); << 142
>> 143 return SandiaCof[0]/energy + SandiaCof[1]/energy2 +
>> 144 SandiaCof[2]/energy3 + SandiaCof[3]/energy4;
140 } 145 }
141 146
142 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
143 148
144 void 149 void
145 G4PEEffectFluoModel::SampleSecondaries(std::ve 150 G4PEEffectFluoModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
146 const G4MaterialCutsCouple* cou 151 const G4MaterialCutsCouple* couple,
147 const G4DynamicParticle* aDynam 152 const G4DynamicParticle* aDynamicPhoton,
148 G4double, 153 G4double,
149 G4double) 154 G4double)
150 { 155 {
151 SetCurrentCouple(couple); 156 SetCurrentCouple(couple);
152 const G4Material* aMaterial = couple->GetMat 157 const G4Material* aMaterial = couple->GetMaterial();
153 158
154 G4double energy = aDynamicPhoton->GetKinetic 159 G4double energy = aDynamicPhoton->GetKineticEnergy();
155 160
156 // select randomly one element constituing t 161 // select randomly one element constituing the material.
157 const G4Element* anElement = SelectRandomAto 162 const G4Element* anElement = SelectRandomAtom(aMaterial,theGamma,energy);
158 163
159 // 164 //
160 // Photo electron 165 // Photo electron
161 // 166 //
162 167
163 // Select atomic shell 168 // Select atomic shell
164 G4int nShells = anElement->GetNbOfAtomicShel 169 G4int nShells = anElement->GetNbOfAtomicShells();
165 G4int i = 0; 170 G4int i = 0;
166 for(; i<nShells; ++i) { 171 for(; i<nShells; ++i) {
167 /* 172 /*
168 G4cout << "i= " << i << " E(eV)= " << ener 173 G4cout << "i= " << i << " E(eV)= " << energy/eV
169 << " Eb(eV)= " << anElement->GetAtomi 174 << " Eb(eV)= " << anElement->GetAtomicShell(i)/eV
170 << " " << anElement->GetName() 175 << " " << anElement->GetName()
171 << G4endl; 176 << G4endl;
172 */ 177 */
173 if(energy >= anElement->GetAtomicShell(i)) 178 if(energy >= anElement->GetAtomicShell(i)) { break; }
174 } 179 }
175 180
176 G4double edep = energy; 181 G4double edep = energy;
177 182
178 // photo-electron is not sampled if shell is << 183 // Normally one shell is available
179 // the flag of photoeffect is "false" and sh << 184 if (i < nShells) {
180 if ( (fPEBelowKShell || 0 == i) && i < nShel << 185
181 <<
182 G4double bindingEnergy = anElement->GetAto 186 G4double bindingEnergy = anElement->GetAtomicShell(i);
183 edep = bindingEnergy; 187 edep = bindingEnergy;
184 G4double esec = 0.0; 188 G4double esec = 0.0;
185 189
186 // sample deexcitation cascade << 190 // sample deexcitation
187 // 191 //
188 if(nullptr != fAtomDeexcitation) { << 192 if(fAtomDeexcitation) {
189 G4int index = couple->GetIndex(); 193 G4int index = couple->GetIndex();
190 if(fAtomDeexcitation->CheckDeexcitationA 194 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
191 G4int Z = G4lrint(anElement->GetZ()); 195 G4int Z = G4lrint(anElement->GetZ());
192 auto as = (G4AtomicShellEnumerator)(i); << 196 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
193 const G4AtomicShell* shell = fAtomDeexcitati 197 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
194 G4double eshell = shell->BindingEnergy 198 G4double eshell = shell->BindingEnergy();
195 if(eshell > bindingEnergy && eshell <= 199 if(eshell > bindingEnergy && eshell <= energy) {
196 bindingEnergy = eshell; 200 bindingEnergy = eshell;
197 edep = eshell; 201 edep = eshell;
198 } 202 }
199 std::size_t nbefore = fvect->size(); << 203 G4int nbefore = fvect->size();
200 fAtomDeexcitation->GenerateParticles(fvect, 204 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
201 std::size_t nafter = fvect->size(); << 205 G4int nafter = fvect->size();
202 for (std::size_t j=nbefore; j<nafter; ++j) { << 206 for (G4int j=nbefore; j<nafter; ++j) {
203 G4double e = ((*fvect)[j])->GetKineticEner 207 G4double e = ((*fvect)[j])->GetKineticEnergy();
204 if(esec + e > edep) { 208 if(esec + e > edep) {
205 // correct energy in order to have energ 209 // correct energy in order to have energy balance
206 e = edep - esec; 210 e = edep - esec;
207 ((*fvect)[j])->SetKineticEnergy(e); 211 ((*fvect)[j])->SetKineticEnergy(e);
208 esec += e; 212 esec += e;
209 /* 213 /*
210 G4cout << "### G4PEffectFluoModel Edep 214 G4cout << "### G4PEffectFluoModel Edep(eV)= " << edep/eV
211 << " Esec(eV)= " << esec/eV 215 << " Esec(eV)= " << esec/eV
212 << " E["<< j << "](eV)= " << e/eV 216 << " E["<< j << "](eV)= " << e/eV
213 << " N= " << nafter 217 << " N= " << nafter
214 << " Z= " << Z << " shell= " << i 218 << " Z= " << Z << " shell= " << i
215 << " Ebind(keV)= " << bindingEnergy/ 219 << " Ebind(keV)= " << bindingEnergy/keV
216 << " Eshell(keV)= " << shell->Bindin 220 << " Eshell(keV)= " << shell->BindingEnergy()/keV
217 << G4endl; 221 << G4endl;
218 */ 222 */
219 // delete the rest of secondaries (shoul 223 // delete the rest of secondaries (should not happens)
220 for (std::size_t jj=nafter-1; jj>j; --jj << 224 for (G4int jj=nafter-1; jj>j; --jj) {
221 delete (*fvect)[jj]; 225 delete (*fvect)[jj];
222 fvect->pop_back(); 226 fvect->pop_back();
223 } 227 }
224 break; 228 break;
225 } 229 }
226 esec += e; 230 esec += e;
227 } 231 }
228 edep -= esec; 232 edep -= esec;
229 } 233 }
230 } 234 }
231 // create photo electron 235 // create photo electron
232 // 236 //
233 G4double elecKineEnergy = energy - binding 237 G4double elecKineEnergy = energy - bindingEnergy;
234 if (elecKineEnergy > fminimalEnergy) { 238 if (elecKineEnergy > fminimalEnergy) {
235 auto aParticle = new G4DynamicParticle(t << 239 G4DynamicParticle* aParticle = new G4DynamicParticle(theElectron,
236 GetAngularDistribution()->SampleDirection(aD 240 GetAngularDistribution()->SampleDirection(aDynamicPhoton,
237 elecKineEnergy, 241 elecKineEnergy,
238 i, couple->GetMaterial()), 242 i, couple->GetMaterial()),
239 elecKineEnergy); 243 elecKineEnergy);
240 fvect->push_back(aParticle); 244 fvect->push_back(aParticle);
241 } else { 245 } else {
242 edep += elecKineEnergy; 246 edep += elecKineEnergy;
243 elecKineEnergy = 0.0; 247 elecKineEnergy = 0.0;
244 } 248 }
245 if(std::abs(energy - elecKineEnergy - esec 249 if(std::abs(energy - elecKineEnergy - esec - edep) > CLHEP::eV) {
246 G4cout << "### G4PEffectFluoModel dE(eV) 250 G4cout << "### G4PEffectFluoModel dE(eV)= "
247 << (energy - elecKineEnergy - esec - ed 251 << (energy - elecKineEnergy - esec - edep)/eV
248 << " shell= " << i 252 << " shell= " << i
249 << " E(keV)= " << energy/keV 253 << " E(keV)= " << energy/keV
250 << " Ebind(keV)= " << bindingEnergy/ke 254 << " Ebind(keV)= " << bindingEnergy/keV
251 << " Ee(keV)= " << elecKineEnergy/keV 255 << " Ee(keV)= " << elecKineEnergy/keV
252 << " Esec(keV)= " << esec/keV 256 << " Esec(keV)= " << esec/keV
253 << " Edep(keV)= " << edep/keV 257 << " Edep(keV)= " << edep/keV
254 << G4endl; 258 << G4endl;
255 } 259 }
256 } 260 }
257 261
258 // kill primary photon 262 // kill primary photon
259 fParticleChange->SetProposedKineticEnergy(0. 263 fParticleChange->SetProposedKineticEnergy(0.);
260 fParticleChange->ProposeTrackStatus(fStopAnd 264 fParticleChange->ProposeTrackStatus(fStopAndKill);
261 if(edep > 0.0) { 265 if(edep > 0.0) {
262 fParticleChange->ProposeLocalEnergyDeposit 266 fParticleChange->ProposeLocalEnergyDeposit(edep);
263 } 267 }
264 } 268 }
265 269
266 //....oooOO0OOooo........oooOO0OOooo........oo 270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
267 271