Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PEEffectFluoModel.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4PEEffectFluoModel.cc 93362 2015-10-19 13:45:19Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4PEEffectFluoModel
34 //
35 // Author: Vladimir Ivanchenko on base of G4PEEffectModel
36 //
37 // Creation date: 13.06.2010
38 //
39 // Modifications:
40 //
41 // Class Description:
42 // Implementation of the photo-electric effect with deexcitation
43 //
44 // -------------------------------------------------------------------
45 //
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 
49 #include "G4PEEffectFluoModel.hh"
50 #include "G4PhysicalConstants.hh"
51 #include "G4SystemOfUnits.hh"
52 #include "G4Electron.hh"
53 #include "G4Gamma.hh"
54 #include "Randomize.hh"
55 #include "G4Material.hh"
56 #include "G4DataVector.hh"
58 #include "G4VAtomDeexcitation.hh"
59 #include "G4LossTableManager.hh"
61 
62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 
64 using namespace std;
65 
67  : G4VEmModel(nam)
68 {
69  theGamma = G4Gamma::Gamma();
70  theElectron = G4Electron::Electron();
71  fminimalEnergy = 1.0*eV;
72  SetDeexcitationFlag(true);
73  fParticleChange = nullptr;
74  fAtomDeexcitation = nullptr;
75 
76  fSandiaCof.resize(4,0.0);
77 
78  // default generator
80 }
81 
82 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
83 
85 {}
86 
87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
88 
90  const G4DataVector&)
91 {
92  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
93  if(nullptr == fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
94  size_t nmat = G4Material::GetNumberOfMaterials();
95  fMatEnergyTh.resize(nmat, 0.0);
96  for(size_t i=0; i<nmat; ++i) {
97  fMatEnergyTh[i] = (*(G4Material::GetMaterialTable()))[i]
98  ->GetSandiaTable()->GetSandiaCofForMaterial(0, 0);
99  //G4cout << "G4PEEffectFluoModel::Initialise Eth(eV)= "
100  // << fMatEnergyTh[i]/eV << G4endl;
101  }
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
105 
106 G4double
111 {
112  // This method may be used only if G4MaterialCutsCouple pointer
113  // has been set properly
115  ->GetSandiaTable()->GetSandiaCofPerAtom((G4int)Z, energy, fSandiaCof);
116 
117  G4double energy2 = energy*energy;
118  G4double energy3 = energy*energy2;
119  G4double energy4 = energy2*energy2;
120 
121  return fSandiaCof[0]/energy + fSandiaCof[1]/energy2 +
122  fSandiaCof[2]/energy3 + fSandiaCof[3]/energy4;
123 }
124 
125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126 
127 G4double
129  const G4ParticleDefinition*,
132 {
133  // This method may be used only if G4MaterialCutsCouple pointer
134  // has been set properly
135  energy = std::max(energy, fMatEnergyTh[material->GetIndex()]);
136  const G4double* SandiaCof =
137  material->GetSandiaTable()->GetSandiaCofForMaterial(energy);
138 
139  G4double energy2 = energy*energy;
140  G4double energy3 = energy*energy2;
141  G4double energy4 = energy2*energy2;
142 
143  return SandiaCof[0]/energy + SandiaCof[1]/energy2 +
144  SandiaCof[2]/energy3 + SandiaCof[3]/energy4;
145 }
146 
147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
148 
149 void
150 G4PEEffectFluoModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
151  const G4MaterialCutsCouple* couple,
152  const G4DynamicParticle* aDynamicPhoton,
153  G4double,
154  G4double)
155 {
156  SetCurrentCouple(couple);
157  const G4Material* aMaterial = couple->GetMaterial();
158 
159  G4double energy = aDynamicPhoton->GetKineticEnergy();
160 
161  // select randomly one element constituing the material.
162  const G4Element* anElement = SelectRandomAtom(aMaterial,theGamma,energy);
163 
164  //
165  // Photo electron
166  //
167 
168  // Select atomic shell
169  G4int nShells = anElement->GetNbOfAtomicShells();
170  G4int i = 0;
171  for(; i<nShells; ++i) {
172  /*
173  G4cout << "i= " << i << " E(eV)= " << energy/eV
174  << " Eb(eV)= " << anElement->GetAtomicShell(i)/eV
175  << " " << anElement->GetName()
176  << G4endl;
177  */
178  if(energy >= anElement->GetAtomicShell(i)) { break; }
179  }
180 
181  G4double edep = energy;
182 
183  // Normally one shell is available
184  if (i < nShells) {
185 
186  G4double bindingEnergy = anElement->GetAtomicShell(i);
187  edep = bindingEnergy;
188  G4double esec = 0.0;
189 
190  // sample deexcitation
191  //
192  if(fAtomDeexcitation) {
193  G4int index = couple->GetIndex();
194  if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
195  G4int Z = G4lrint(anElement->GetZ());
197  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
198  G4double eshell = shell->BindingEnergy();
199  if(eshell > bindingEnergy && eshell <= energy) {
200  bindingEnergy = eshell;
201  edep = eshell;
202  }
203  G4int nbefore = fvect->size();
204  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
205  G4int nafter = fvect->size();
206  for (G4int j=nbefore; j<nafter; ++j) {
207  G4double e = ((*fvect)[j])->GetKineticEnergy();
208  if(esec + e > edep) {
209  // correct energy in order to have energy balance
210  e = edep - esec;
211  ((*fvect)[j])->SetKineticEnergy(e);
212  esec += e;
213  /*
214  G4cout << "### G4PEffectFluoModel Edep(eV)= " << edep/eV
215  << " Esec(eV)= " << esec/eV
216  << " E["<< j << "](eV)= " << e/eV
217  << " N= " << nafter
218  << " Z= " << Z << " shell= " << i
219  << " Ebind(keV)= " << bindingEnergy/keV
220  << " Eshell(keV)= " << shell->BindingEnergy()/keV
221  << G4endl;
222  */
223  // delete the rest of secondaries (should not happens)
224  for (G4int jj=nafter-1; jj>j; --jj) {
225  delete (*fvect)[jj];
226  fvect->pop_back();
227  }
228  break;
229  }
230  esec += e;
231  }
232  edep -= esec;
233  }
234  }
235  // create photo electron
236  //
237  G4double elecKineEnergy = energy - bindingEnergy;
238  if (elecKineEnergy > fminimalEnergy) {
239  G4DynamicParticle* aParticle = new G4DynamicParticle(theElectron,
240  GetAngularDistribution()->SampleDirection(aDynamicPhoton,
241  elecKineEnergy,
242  i, couple->GetMaterial()),
243  elecKineEnergy);
244  fvect->push_back(aParticle);
245  } else {
246  edep += elecKineEnergy;
247  elecKineEnergy = 0.0;
248  }
249  if(fabs(energy - elecKineEnergy - esec - edep) > eV) {
250  G4cout << "### G4PEffectFluoModel dE(eV)= "
251  << (energy - elecKineEnergy - esec - edep)/eV
252  << " shell= " << i
253  << " E(keV)= " << energy/keV
254  << " Ebind(keV)= " << bindingEnergy/keV
255  << " Ee(keV)= " << elecKineEnergy/keV
256  << " Esec(keV)= " << esec/keV
257  << " Edep(keV)= " << edep/keV
258  << G4endl;
259  }
260  }
261 
262  // kill primary photon
263  fParticleChange->SetProposedKineticEnergy(0.);
264  fParticleChange->ProposeTrackStatus(fStopAndKill);
265  if(edep > 0.0) {
266  fParticleChange->ProposeLocalEnergyDeposit(edep);
267  }
268 }
269 
270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4int GetNbOfAtomicShells() const
Definition: G4Element.hh:147
G4double GetKineticEnergy() const
size_t GetIndex() const
Definition: G4Material.hh:262
G4double GetSandiaCofForMaterial(G4int, G4int) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) override
G4double GetZ() const
Definition: G4Element.hh:131
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:587
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:619
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4SandiaTable * GetSandiaTable() const
Definition: G4Material.hh:229
string material
Definition: eplot.py:19
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
const G4MaterialCutsCouple * CurrentCouple() const
Definition: G4VEmModel.hh:454
G4GLOB_DLL std::ostream G4cout
static constexpr double eV
Definition: G4SIunits.hh:215
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy) override
static size_t GetNumberOfMaterials()
Definition: G4Material.cc:594
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A, G4double, G4double) override
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:447
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:626
static G4Electron * Electron()
Definition: G4Electron.cc:94
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
G4double GetAtomicShell(G4int index) const
Definition: G4Element.cc:367
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788
G4double bindingEnergy(G4int A, G4int Z)
static constexpr double keV
Definition: G4SIunits.hh:216
G4PEEffectFluoModel(const G4String &nam="PhotoElectric")
void GetSandiaCofPerAtom(G4int Z, G4double energy, std::vector< G4double > &coeff) const
G4AtomicShellEnumerator
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:132
const G4Material * GetMaterial() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &) override