source: trunk/source/processes/electromagnetic/lowenergy/src/G4ShellData.cc@ 1197

Last change on this file since 1197 was 1196, checked in by garnier, 16 years ago

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