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