ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4RDShellData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4RDShellData.cc
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 
53  : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
54 { }
55 
56 // Destructor
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 
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 
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 
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 
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 
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  char* path = std::getenv("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 
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 
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 }