Geant4  10.02.p01
G4RDShellData.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$
28 // GEANT4 tag $Name: $
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 "G4RDShellData.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 
55  : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
56 { }
57 
58 // Destructor
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 
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 
100 const std::vector<G4double>& G4RDShellData::ShellIdVector(G4int Z) const
101 {
102  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
103  if (Z < zMin || Z > zMax)
104  G4Exception("G4RDShellData::ShellIdVector()", "OutOfRange",
105  FatalException, "Z outside boundaries!");
106  pos = idMap.find(Z);
107  std::vector<G4double>* dataSet = (*pos).second;
108  return *dataSet;
109 }
110 
111 
112 const std::vector<G4double>& G4RDShellData::ShellVector(G4int Z) const
113 {
114  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
115  if (Z < zMin || Z > zMax)
116  G4Exception("G4RDShellData::ShellVector()", "OutOfRange",
117  FatalException, "Z outside boundaries!");
118  pos = occupancyPdfMap.find(Z);
119  std::vector<G4double>* dataSet = (*pos).second;
120  return *dataSet;
121 }
122 
123 
125 {
126  G4int n = -1;
127 
128  if (Z >= zMin && Z <= zMax)
129  {
130  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
131  pos = idMap.find(Z);
132  if (pos!= idMap.end())
133  {
134  std::vector<G4double> dataSet = *((*pos).second);
135  G4int nData = dataSet.size();
136  if (shellIndex >= 0 && shellIndex < nData)
137  {
138  n = (G4int) dataSet[shellIndex];
139  }
140  }
141  }
142  return n;
143 }
144 
145 
147 {
148  G4double prob = -1.;
149 
150  if (Z >= zMin && Z <= zMax)
151  {
152  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
153  pos = idMap.find(Z);
154  if (pos!= idMap.end())
155  {
156  std::vector<G4double> dataSet = *((*pos).second);
157  G4int nData = dataSet.size();
158  if (shellIndex >= 0 && shellIndex < nData)
159  {
160  prob = dataSet[shellIndex];
161  }
162  }
163  }
164  return prob;
165 }
166 
167 
168 
170 {
171  G4double value = 0.;
172 
173  if (Z >= zMin && Z <= zMax)
174  {
175  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
176  pos = bindingMap.find(Z);
177  if (pos!= bindingMap.end())
178  {
179  G4DataVector dataSet = *((*pos).second);
180  G4int nData = dataSet.size();
181  if (shellIndex >= 0 && shellIndex < nData)
182  {
183  value = dataSet[shellIndex];
184  }
185  }
186  }
187  return value;
188 }
189 
191 {
192  for (G4int Z = zMin; Z <= zMax; Z++)
193  {
194  G4cout << "---- Shell data for Z = "
195  << Z
196  << " ---- "
197  << G4endl;
198  G4int nSh = nShells[Z-1];
199  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posId;
200  posId = idMap.find(Z);
201  std::vector<G4double>* ids = (*posId).second;
202  std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posE;
203  posE = bindingMap.find(Z);
204  G4DataVector* energies = (*posE).second;
205  for (G4int i=0; i<nSh; i++)
206  {
207  G4int id = (G4int) (*ids)[i];
208  G4double e = (*energies)[i] / keV;
209  G4cout << i << ") ";
210 
211  if (occupancyData)
212  {
213  G4cout << " Occupancy: ";
214  }
215  else
216  {
217  G4cout << " Shell id: ";
218  }
219  G4cout << id << " - Binding energy = "
220  << e << " keV ";
221  if (occupancyData)
222  {
223  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posOcc;
224  posOcc = occupancyPdfMap.find(Z);
225  std::vector<G4double> probs = *((*posOcc).second);
226  G4double prob = probs[i];
227  G4cout << "- Probability = " << prob;
228  }
229  G4cout << G4endl;
230  }
231  G4cout << "-------------------------------------------------"
232  << G4endl;
233  }
234 }
235 
236 
237 void G4RDShellData::LoadData(const G4String& fileName)
238 {
239  // Build the complete string identifying the file with the data set
240 
241  std::ostringstream ost;
242 
243  ost << fileName << ".dat";
244 
245  G4String name(ost.str());
246 
247  char* path = getenv("G4LEDATA");
248  if (!path)
249  {
250  G4String excep("G4LEDATA environment variable not set!");
251  G4Exception("G4RDShellData::LoadData()", "InvalidSetup",
252  FatalException, excep);
253  }
254 
255  G4String pathString(path);
256  G4String dirFile = pathString + name;
257  std::ifstream file(dirFile);
258  std::filebuf* lsdp = file.rdbuf();
259 
260  if (! (lsdp->is_open()) )
261  {
262  G4String s1("Data file: ");
263  G4String s2(" not found");
264  G4String excep = s1 + dirFile + s2;
265  G4Exception("G4RDShellData::LoadData()", "DataNotFound",
266  FatalException, excep);
267  }
268 
269  G4double a = 0;
270  G4int k = 1;
271  G4int s = 0;
272 
273  G4int Z = 1;
274  G4DataVector* energies = new G4DataVector;
275  std::vector<G4double>* ids = new std::vector<G4double>;
276 
277  do {
278  file >> a;
279  G4int nColumns = 2;
280  if (a == -1)
281  {
282  if (s == 0)
283  {
284  // End of a shell data set
285  idMap[Z] = ids;
286  bindingMap[Z] = energies;
287  G4int n = ids->size();
288  nShells.push_back(n);
289  // Start of new shell data set
290  ids = new std::vector<G4double>;
291  energies = new G4DataVector;
292  Z++;
293  }
294  s++;
295  if (s == nColumns)
296  {
297  s = 0;
298  }
299  }
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 
326  // For Doppler broadening: the data set contains shell occupancy and binding energy for each shell
327  // Build additional map with probability for each shell based on its occupancy
328 
329  if (occupancyData)
330  {
331  // Build cumulative from raw shell occupancy
332 
333  for (G4int Z=zMin; Z <= zMax; Z++)
334  {
335  std::vector<G4double> occupancy = ShellIdVector(Z);
336 
337  std::vector<G4double>* prob = new std::vector<G4double>;
338  G4double scale = 1. / G4double(Z);
339 
340  prob->push_back(occupancy[0] * scale);
341  for (size_t i=1; i<occupancy.size(); i++)
342  {
343  prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
344  }
345  occupancyPdfMap[Z] = prob;
346 
347  /*
348  G4double scale = 1. / G4double(Z);
349  // transform((*prob).begin(),(*prob).end(),(*prob).begin(),bind2nd(multiplies<G4double>(),scale));
350 
351  for (size_t i=0; i<occupancy.size(); i++)
352  {
353  (*prob)[i] *= scale;
354  }
355  */
356  }
357  }
358 }
359 
360 
362 {
363  if (Z < zMin || Z > zMax)
364  G4Exception("G4RDShellData::SelectRandomShell()", "OutOfRange",
365  FatalException, "Z outside boundaries!");
366 
367  G4int shellIndex = 0;
368  std::vector<G4double> prob = ShellVector(Z);
369  G4double random = G4UniformRand();
370 
371  // std::vector<G4double>::const_iterator pos;
372  // pos = lower_bound(prob.begin(),prob.end(),random);
373 
374  // Binary search the shell with probability less or equal random
375 
377  G4int upperBound = nShells;
378 
379  while (shellIndex <= upperBound)
380  {
381  G4int midShell = (shellIndex + upperBound) / 2;
382  if ( random < prob[midShell] )
383  upperBound = midShell - 1;
384  else
385  shellIndex = midShell + 1;
386  }
387  if (shellIndex >= nShells) shellIndex = nShells - 1;
388 
389  return shellIndex;
390 }
static const double MeV
Definition: G4SIunits.hh:211
std::map< G4int, std::vector< G4double > *, std::less< G4int > > idMap
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
G4double BindingEnergy(G4int Z, G4int shellIndex) const
G4double a
Definition: TRTMaterials.hh:39
int G4int
Definition: G4Types.hh:78
size_t NumberOfShells(G4int Z) const
void LoadData(const G4String &fileName)
static const double s
Definition: G4SIunits.hh:168
void PrintData() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
std::vector< G4int > nShells
bool G4bool
Definition: G4Types.hh:79
const std::vector< G4double > & ShellVector(G4int Z) const
G4RDShellData(G4int minZ=1, G4int maxZ=100, G4bool isOccupancy=false)
const G4int n
G4int ShellId(G4int Z, G4int shellIndex) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const std::vector< G4double > & ShellIdVector(G4int Z) const
std::map< G4int, std::vector< G4double > *, std::less< G4int > > occupancyPdfMap
G4double ShellOccupancyProbability(G4int Z, G4int shellIndex) const
#define G4endl
Definition: G4ios.hh:61
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
G4bool occupancyData
G4int SelectRandomShell(G4int Z) const
static const G4double pos
std::map< G4int, G4DataVector *, std::less< G4int > > bindingMap