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