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