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