Geant4_10
G4PenelopeIonisationCrossSection.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: G4PenelopeIonisationCrossSection.cc 74626 2013-10-17 07:00:59Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 14 Mar 2012 L Pandola First complete implementation for e-
33 //
34 //
37 //
39 #include "G4SystemOfUnits.hh"
43 #include "G4Electron.hh"
45 
47  G4VhShellCrossSection("Penelope"),shellIDTable(0),
48  theCrossSectionHandler(0)
49 {
51  nMaxLevels = 9;
52 
53  // Verbosity scale:
54  // 0 = nothing
55  // 1 = calculation of cross sections, file openings, sampling of atoms
56  // 2 = entering in methods
57  verboseLevel = 0;
58 
59  fLowEnergyLimit = 10.0*eV;
60  fHighEnergyLimit = 100.0*GeV;
61 
62  transitionManager = G4AtomicTransitionManager::Instance();
63 }
64 
65 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
66 
68 {
69  if (theCrossSectionHandler)
70  delete theCrossSectionHandler;
71 }
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
77  G4double incidentEnergy,
78  G4double ,
79  const G4Material* material)
80 {
81  if (verboseLevel > 1)
82  G4cout << "Entering in method G4PenelopeIonisationCrossSection::CrossSection()" << G4endl;
83 
84  G4double cross = 0.;
85 
86  //Material pointer is not available
87  if (!material)
88  {
89  //CRASH!
91  ed << "The method has been called with a NULL G4Material pointer" << G4endl;
92  G4Exception("G4PenelopeIonisationCrossSection::CrossSection()","em2042",
93  FatalException,ed);
94  return cross;
95  }
96 
97  if (!theCrossSectionHandler)
98  theCrossSectionHandler = new G4PenelopeIonisationXSHandler();
99 
100  theCrossSectionHandler->BuildXSTable(material,0.,G4Electron::Electron());
101 
102  G4int nmax = std::min(nMaxLevels,transitionManager->NumberOfShells(Z));
103 
104  if(G4int(shell) < nmax &&
105  incidentEnergy >= fLowEnergyLimit && incidentEnergy <= fHighEnergyLimit)
106  {
107  //The shells in Penelope are organized per *material*, rather than per
108  //element, so given a material one has to find the proper index for the
109  //given Z and shellID. An appropriate lookup table is used to avoid
110  //recalculation.
111  G4int index = FindShellIDIndex(material,Z,shell);
112 
113  //Index is not available!
114  if (index < 0)
115  return cross;
116 
117  const G4PenelopeCrossSection* theXS =
118  theCrossSectionHandler->GetCrossSectionTableForCouple(G4Electron::Electron(),
119  material,
120  0.);
121 
122  //Cross check that everything is fine:
123  G4PenelopeOscillator* theOsc = oscManager->GetOscillatorIonisation(material,index);
124  if (theOsc->GetParentZ() != Z || theOsc->GetShellFlag()-1 != G4int(shell))
125  {
126  //something went wrong!
128  ed << "There is something wrong here: it looks like the index is wrong" << G4endl;
129  ed << "Requested: shell " << G4int(shell) << " and Z = " << Z << G4endl;
130  ed << "Retrieved: " << theOsc->GetShellFlag()-1 << " and Z = " << theOsc->GetParentZ() << G4endl;
131  G4Exception("G4PenelopeIonisationCrossSection::CrossSection()","em2043",
132  JustWarning,ed);
133  return cross;
134  }
135 
136 
137  G4double crossPerMolecule = (theXS) ? theXS->GetShellCrossSection(index,incidentEnergy) : 0.;
138 
139  //Now it must be converted in cross section per atom. I need the number of
140  //atoms of the given Z per molecule.
141  G4double atomsPerMolec = oscManager->GetNumberOfZAtomsPerMolecule(material,Z);
142  if (atomsPerMolec)
143  cross = crossPerMolecule/atomsPerMolec;
144 
145 
146  if (verboseLevel > 0)
147  {
148  G4cout << "Cross section of shell " << G4int(shell) << " and Z= " << Z;
149  G4cout << " of material: " << material->GetName() << " and energy = " << incidentEnergy/keV << " keV" << G4endl;
150  G4cout << "--> " << cross/barn << " barn" << G4endl;
151  G4cout << "Shell binding energy: " << theOsc->GetIonisationEnergy()/eV << " eV;" ;
152  G4cout << " resonance energy: " << theOsc->GetResonanceEnergy()/eV << "eV" << G4endl;
153  if (verboseLevel > 2)
154  {
155  G4cout << "Cross section per molecule: " << crossPerMolecule/barn << " barn" << G4endl;
156  G4cout << "Atoms " << Z << " per molecule: " << atomsPerMolec << G4endl;
157  }
158  }
159  }
160 
161  return cross;
162 }
163 
164 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
165 std::vector<G4double>
167  G4double kinEnergy,
168  G4double, G4double,
169  const G4Material* mat)
170 {
171  G4int nmax = std::min(nMaxLevels,transitionManager->NumberOfShells(Z));
172  std::vector<G4double> vec(nmax,0.0);
173  for(G4int i=0; i<nmax; ++i) {
174  vec[i] = CrossSection(Z, G4AtomicShellEnumerator(i), kinEnergy,0.,mat);
175  }
176  return vec;
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180 
181 std::vector<G4double>
183  G4double kinEnergy,
184  G4double,
185  G4double,
186  const G4Material* mat)
187 {
188  std::vector<G4double> vec = GetCrossSection(Z, kinEnergy,0,0,mat);
189  size_t n = vec.size();
190  size_t i=0;
191  G4double sum = 0.0;
192  for(i=0; i<n; ++i) { sum += vec[i]; }
193  if(sum > 0.0) {
194  sum = 1.0/sum;
195  for(i=0; i<n; ++i) { vec[i] = vec[i]*sum; }
196  }
197  return vec;
198 }
199 
200 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
201 
202 G4int G4PenelopeIonisationCrossSection::FindShellIDIndex(const G4Material* mat,
203  G4int Z,
205 {
206  if (verboseLevel > 1)
207  G4cout << "Entering in method G4PenelopeIonisationCrossSection::FindShellIDIndex()" << G4endl;
208 
209  if (!shellIDTable)
210  shellIDTable = new std::map< std::pair<const G4Material*,G4int>, G4DataVector*>;
211 
212  std::pair<const G4Material*,G4int> theKey = std::make_pair(mat,Z);
213  G4int result = -1;
214  G4int ishell = G4int(shell);
215 
216  if (shellIDTable->count(theKey)) //table already built, and containing the element
217  {
218  if (verboseLevel > 2)
219  G4cout << "FindShellIDIndex: Table already built for " << mat->GetName() << G4endl;
220  G4DataVector* theVec = shellIDTable->find(theKey)->second;
221 
222  if (ishell>=0 && ishell < (G4int) theVec->size()) //check we are not off-boundary
223  result = (G4int) (*theVec)[ishell];
224  else
225  {
227  ed << "Shell ID: " << ishell << " not available for material " << mat->GetName() << " and Z = " <<
228  Z << G4endl;
229  G4Exception("G4PenelopeIonisationCrossSection::FindShellIDIndex()","em2041",JustWarning,
230  ed);
231  return -1;
232  }
233  }
234  else
235  {
236  if (verboseLevel > 2)
237  G4cout << "FindShellIDIndex: Table to be built for " << mat->GetName() << G4endl;
238  //Not contained: look for it
239  G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(mat);
240  size_t numberOfOscillators = theTable->size();
241 
242  //oscillator loop
243  //initialize everything at -1
244  G4DataVector* dat = new G4DataVector(nMaxLevels,-1);
245  for (size_t iosc=0;iosc<numberOfOscillators;iosc++)
246  {
247  G4PenelopeOscillator* theOsc = (*theTable)[iosc];
248  //level is found!
249  if (theOsc->GetParentZ() == Z)
250  {
251  //individual shells relative to the given material
252  G4int shFlag = theOsc->GetShellFlag();
253  //Notice: GetShellFlag() starts from 1, the G4AtomicShellEnumerator from 0
254  if (shFlag < 30)
255  (*dat)[shFlag-1] = (G4double) iosc; //index of the given shell
256  if ((shFlag-1) == ishell) // this is what we were looking for
257  result = (G4int) iosc;
258  }
259  }
260  shellIDTable->insert(std::make_pair(theKey,dat));
261  }
262 
263  if (verboseLevel > 1)
264  G4cout << "Leaving method G4PenelopeIonisationCrossSection::FindShellIDIndex() with index = " << result << G4endl;
265 
266  return result;
267 
268 }
std::vector< G4double > Probabilities(G4int Z, G4double incidentEnergy, G4double mass, G4double deltaEnergy, const G4Material *mat)
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
Int_t index
Definition: macro.C:9
const G4String & GetName() const
Definition: G4Material.hh:176
G4double G4NeutronHPJENDLHEData::G4double result
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
G4PenelopeOscillator * GetOscillatorIonisation(const G4Material *, G4int)
int G4int
Definition: G4Types.hh:78
Float_t mat
Definition: plot.C:40
string material
Definition: eplot.py:19
Char_t n[5]
G4GLOB_DLL std::ostream G4cout
Float_t Z
Definition: plot.C:39
const G4int nmax
static G4PenelopeOscillatorManager * GetOscillatorManager()
G4double GetResonanceEnergy() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
~G4PenelopeIonisationCrossSection()
Destructor. Clean all tables.
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
std::vector< G4double > GetCrossSection(G4int Z, G4double incidentEnergy, G4double mass, G4double deltaEnergy, const G4Material *mat)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
double G4double
Definition: G4Types.hh:76
G4double CrossSection(G4int Z, G4AtomicShellEnumerator shell, G4double incidentEnergy, G4double mass, const G4Material *mat)
G4double GetNumberOfZAtomsPerMolecule(const G4Material *, G4int Z)
static G4AtomicTransitionManager * Instance()
G4AtomicShellEnumerator
G4double GetShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (per molecule)