| 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 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // -----------
32 // 31 Jul 2001 MGP Created
33 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings
34 //
35 // -------------------------------------------------------------------
36
37 #include "G4ShellData.hh"
38 #include "G4DataVector.hh"
39 #include "G4SystemOfUnits.hh"
40 #include <fstream>
41 #include <sstream>
42 #include <numeric>
43 #include <algorithm>
44 #include <valarray>
45 #include <functional>
46 #include "Randomize.hh"
47
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49
50 G4ShellData::G4ShellData(G4int minZ, G4int maxZ, G4bool isOccupancy)
51 : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
52 {}
53
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56 G4ShellData::~G4ShellData()
57 {
58 for (auto& pos : idMap)
59 {
60 std::vector<G4double>* dataSet = pos.second;
61 delete dataSet;
62 }
63 for (auto& pos2 : bindingMap)
64 {
65 G4DataVector* dataSet = pos2.second;
66 delete dataSet;
67 }
68
69 if (occupancyData)
70 {
71 for (auto& pos3 : occupancyPdfMap)
72 {
73 std::vector<G4double>* dataSet = pos3.second;
74 delete dataSet;
75 }
76 }
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81 std::size_t G4ShellData::NumberOfShells(G4int Z) const
82 {
83 G4int z = Z - 1;
84 G4int n = 0;
85
86 if (Z>= zMin && Z <= zMax)
87 {
88 n = nShells[z];
89 }
90 return n;
91 }
92
93 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
94
95 const std::vector<G4double>& G4ShellData::ShellIdVector(G4int Z) const
96 {
97 if (Z < zMin || Z > zMax) {
98 G4Exception("G4ShellData::ShellIdVector","de0001",FatalErrorInArgument, "Z outside boundaries");
99 }
100 auto pos = idMap.find(Z);
101 std::vector<G4double>* dataSet = (*pos).second;
102 return *dataSet;
103 }
104
105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106
107 const std::vector<G4double>& G4ShellData::ShellVector(G4int Z) const
108 {
109 if (Z < zMin || Z > zMax)
110 G4Exception("G4ShellData::ShellVector()","de0001",JustWarning,"Z outside boundaries");
111 auto pos = occupancyPdfMap.find(Z);
112 std::vector<G4double>* dataSet = (*pos).second;
113 return *dataSet;
114 }
115
116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117
118 G4int G4ShellData::ShellId(G4int Z, G4int shellIndex) const
119 {
120 G4int n = -1;
121
122 if (Z >= zMin && Z <= zMax)
123 {
124 auto pos = idMap.find(Z);
125 if (pos!= idMap.end())
126 {
127 std::vector<G4double> dataSet = *((*pos).second);
128 G4int nData = (G4int)dataSet.size();
129 if (shellIndex >= 0 && shellIndex < nData)
130 {
131 n = (G4int) dataSet[shellIndex];
132 }
133 }
134 }
135 return n;
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139
140 G4double G4ShellData::ShellOccupancyProbability(G4int Z, G4int shellIndex) const
141 {
142 G4double prob = -1.;
143 if (Z >= zMin && Z <= zMax)
144 {
145 auto pos = idMap.find(Z);
146 if (pos!= idMap.end())
147 {
148 std::vector<G4double> dataSet = *((*pos).second);
149 G4int nData = (G4int)dataSet.size();
150 if (shellIndex >= 0 && shellIndex < nData)
151 {
152 prob = dataSet[shellIndex];
153 }
154 }
155 }
156 return prob;
157 }
158
159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160
161 G4double G4ShellData::BindingEnergy(G4int Z, G4int shellIndex) const
162 {
163 G4double value = 0.;
164
165 if (Z >= zMin && Z <= zMax)
166 {
167 auto pos = bindingMap.find(Z);
168 if (pos!= bindingMap.end())
169 {
170 G4DataVector dataSet = *((*pos).second);
171 G4int nData = (G4int)dataSet.size();
172 if (shellIndex >= 0 && shellIndex < nData)
173 {
174 value = dataSet[shellIndex];
175 }
176 }
177 }
178 return value;
179 }
180
181 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
182
183 void G4ShellData::PrintData() const
184 {
185 for (G4int Z = zMin; Z <= zMax; Z++)
186 {
187 G4cout << "---- Shell data for Z = "
188 << Z
189 << " ---- "
190 << G4endl;
191 G4int nSh = nShells[Z-1];
192 auto posId = idMap.find(Z);
193 std::vector<G4double>* ids = (*posId).second;
194 auto posE = bindingMap.find(Z);
195 G4DataVector* energies = (*posE).second;
196 for (G4int i=0; i<nSh; ++i)
197 {
198 G4int id = (G4int) (*ids)[i];
199 G4double e = (*energies)[i] / keV;
200 G4cout << i << ") ";
201
202 if (occupancyData)
203 {
204 G4cout << " Occupancy: ";
205 }
206 else
207 {
208 G4cout << " Shell id: ";
209 }
210 G4cout << id << " - Binding energy = "
211 << e << " keV ";
212 if (occupancyData)
213 {
214 auto posOcc = occupancyPdfMap.find(Z);
215 std::vector<G4double> probs = *((*posOcc).second);
216 G4double prob = probs[i];
217 G4cout << "- Probability = " << prob;
218 }
219 G4cout << G4endl;
220 }
221 G4cout << "-------------------------------------------------"
222 << G4endl;
223 }
224 }
225
226 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
227
228 void G4ShellData::LoadData(const G4String& fileName)
229 {
230 // Build the complete string identifying the file with the data set
231 std::ostringstream ost;
232 ost << fileName << ".dat";
233 G4String name(ost.str());
234
235 const char* path = G4FindDataDir("G4LEDATA");
236 if (!path)
237 {
238 G4String excep("G4ShellData::LoadData()");
239 G4Exception(excep,"em0006",FatalException,"Please set G4LEDATA");
240 return;
241 }
242
243 G4String pathString(path);
244 G4String dirFile = pathString + name;
245 std::ifstream file(dirFile);
246 std::filebuf* lsdp = file.rdbuf();
247
248 if (! (lsdp->is_open()) )
249 {
250
251 G4String excep = "G4ShellData::LoadData()";
252 G4String msg = "data file: " + dirFile + " not found";
253 G4Exception(excep, "em0003",FatalException, msg );
254 return;
255 }
256
257 G4double a = 0;
258 G4int k = 1;
259 G4int sLocal = 0;
260
261 G4int Z = 1;
262 G4DataVector* energies = new G4DataVector;
263 std::vector<G4double>* ids = new std::vector<G4double>;
264
265 do {
266 file >> a;
267 G4int nColumns = 2;
268 if (a == -1)
269 {
270 if (sLocal == 0)
271 {
272 // End of a shell data set
273 idMap[Z] = ids;
274 bindingMap[Z] = energies;
275 G4int n = (G4int)ids->size();
276 nShells.push_back(n);
277 // Start of new shell data set
278
279 ids = new std::vector<G4double>;
280 energies = new G4DataVector;
281 Z++;
282 }
283 sLocal++;
284 if (sLocal == nColumns)
285 {
286 sLocal = 0;
287 }
288 }
289
290 // moved out of the do-while since might go to a leak.
291 // else if (a == -2)
292 // {
293 // End of file; delete the empty vectors created when encountering the last -1 -1 row
294 // delete energies;
295 // delete ids;
296 //nComponents = components.size();
297 // }
298 else
299 {
300 // 1st column is shell id
301 if(k%nColumns != 0)
302 {
303 ids->push_back(a);
304 k++;
305 }
306 else if (k%nColumns == 0)
307 {
308 // 2nd column is binding energy
309 G4double e = a * MeV;
310 energies->push_back(e);
311 k = 1;
312 }
313 }
314 } while (a != -2); // end of file
315 file.close();
316 delete energies;
317 delete ids;
318
319 // For Doppler broadening: the data set contains shell occupancy and binding energy for each shell
320 // Build additional map with probability for each shell based on its occupancy
321
322 if (occupancyData)
323 {
324 // Build cumulative from raw shell occupancy
325 for (G4int ZLocal=zMin; ZLocal <= zMax; ++ZLocal)
326 {
327 std::vector<G4double> occupancy = ShellIdVector(ZLocal);
328
329 std::vector<G4double>* prob = new std::vector<G4double>;
330 G4double scale = 1. / G4double(ZLocal);
331
332 prob->push_back(occupancy[0] * scale);
333 for (std::size_t i=1; i<occupancy.size(); ++i)
334 {
335 prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
336 }
337 occupancyPdfMap[ZLocal] = prob;
338 }
339 }
340 }
341
342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
343
344 G4int G4ShellData::SelectRandomShell(G4int Z) const
345 {
346 if (Z < zMin || Z > zMax)
347 G4Exception("G4ShellData::SelectrandomShell","de0001",FatalErrorInArgument, "Z outside boundaries");
348
349 G4int shellIndex = 0;
350 std::vector<G4double> prob = ShellVector(Z);
351 G4double random = G4UniformRand();
352
353 // Binary search the shell with probability less or equal random
354 G4int nShellsLocal = (G4int)NumberOfShells(Z);
355 G4int upperBound = nShellsLocal;
356
357 while (shellIndex <= upperBound)
358 {
359 G4int midShell = (shellIndex + upperBound) / 2;
360 if ( random < prob[midShell] )
361 upperBound = midShell - 1;
362 else
363 shellIndex = midShell + 1;
364 }
365 if (shellIndex >= nShellsLocal) shellIndex = nShellsLocal - 1;
366
367 return shellIndex;
368 }
369