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