ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ShellData.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ShellData.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 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings
34 //
35 // -------------------------------------------------------------------
36 
37 #include "G4ShellData.hh"
38 #include "G4DataVector.hh"
39 #include "G4SystemOfUnits.hh"
40 #include <fstream>
41 #include <sstream>
42 #include <numeric>
43 #include <algorithm>
44 #include <valarray>
45 #include <functional>
46 #include "Randomize.hh"
47 
48 // Constructor
49 
51  : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
52 { }
53 
54 // Destructor
56 {
57  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos;
58  for (pos = idMap.begin(); pos != idMap.end(); ++pos)
59  {
60  std::vector<G4double>* dataSet = (*pos).second;
61  delete dataSet;
62  }
63 
64  std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos2;
65  for (pos2 = bindingMap.begin(); pos2 != bindingMap.end(); ++pos2)
66  {
67  G4DataVector* dataSet = (*pos2).second;
68  delete dataSet;
69  }
70 
71  if (occupancyData)
72  {
73  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos3;
74  for (pos3 = occupancyPdfMap.begin(); pos3 != occupancyPdfMap.end(); ++pos3)
75  {
76  std::vector<G4double>* dataSet = (*pos3).second;
77  delete dataSet;
78  }
79  }
80 }
81 
82 
84 {
85  G4int z = Z - 1;
86  G4int n = 0;
87 
88  if (Z>= zMin && Z <= zMax)
89  {
90  n = nShells[z];
91  }
92  return n;
93 }
94 
95 
96 const std::vector<G4double>& G4ShellData::ShellIdVector(G4int Z) const
97 {
98  std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
99  if (Z < zMin || Z > zMax) {
100 
101  G4Exception("G4ShellData::ShellIdVector","de0001",FatalErrorInArgument, "Z outside boundaries");
102 
103 
104  } pos = idMap.find(Z);
105  std::vector<G4double>* dataSet = (*pos).second;
106  return *dataSet;
107 }
108 
109 
110 const 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()","de0001",JustWarning,"Z outside boundaries");
114  pos = occupancyPdfMap.find(Z);
115  std::vector<G4double>* dataSet = (*pos).second;
116  return *dataSet;
117 }
118 
119 
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 
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 
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 
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 
233 void 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 = std::getenv("G4LEDATA");
244  if (!path)
245  {
246  G4String excep("G4ShellData::LoadData()");
247  G4Exception(excep,"em0006",FatalException,"Please set G4LEDATA");
248  return;
249  }
250 
251  G4String pathString(path);
252  G4String dirFile = pathString + name;
253  std::ifstream file(dirFile);
254  std::filebuf* lsdp = file.rdbuf();
255 
256  if (! (lsdp->is_open()) )
257  {
258 
259  G4String excep = "G4ShellData::LoadData()";
260  G4String msg = "data file: " + dirFile + " not found";
261  G4Exception(excep, "em0003",FatalException, msg );
262  return;
263  }
264 
265  G4double a = 0;
266  G4int k = 1;
267  G4int sLocal = 0;
268 
269  G4int Z = 1;
270  G4DataVector* energies = new G4DataVector;
271  std::vector<G4double>* ids = new std::vector<G4double>;
272 
273  do {
274  file >> a;
275  G4int nColumns = 2;
276  if (a == -1)
277  {
278  if (sLocal == 0)
279  {
280  // End of a shell data set
281  idMap[Z] = ids;
282  bindingMap[Z] = energies;
283  G4int n = ids->size();
284  nShells.push_back(n);
285  // Start of new shell data set
286 
287  ids = new std::vector<G4double>;
288  energies = new G4DataVector;
289  Z++;
290  }
291  sLocal++;
292  if (sLocal == nColumns)
293  {
294  sLocal = 0;
295  }
296  }
297 
298  // moved out of the do-while since might go to a leak.
299  // else if (a == -2)
300  // {
301  // End of file; delete the empty vectors created when encountering the last -1 -1 row
302  // delete energies;
303  // delete ids;
304  //nComponents = components.size();
305  // }
306  else
307  {
308  // 1st column is shell id
309  if(k%nColumns != 0)
310  {
311  ids->push_back(a);
312  k++;
313  }
314  else if (k%nColumns == 0)
315  {
316  // 2nd column is binding energy
317  G4double e = a * MeV;
318  energies->push_back(e);
319  k = 1;
320  }
321  }
322  } while (a != -2); // end of file
323  file.close();
324  delete energies;
325  delete ids;
326 
327  // For Doppler broadening: the data set contains shell occupancy and binding energy for each shell
328  // Build additional map with probability for each shell based on its occupancy
329 
330  if (occupancyData)
331  {
332  // Build cumulative from raw shell occupancy
333 
334  for (G4int ZLocal=zMin; ZLocal <= zMax; ZLocal++)
335  {
336  std::vector<G4double> occupancy = ShellIdVector(ZLocal);
337 
338  std::vector<G4double>* prob = new std::vector<G4double>;
339  G4double scale = 1. / G4double(ZLocal);
340 
341  prob->push_back(occupancy[0] * scale);
342  for (size_t i=1; i<occupancy.size(); i++)
343  {
344  prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
345  }
346  occupancyPdfMap[ZLocal] = prob;
347 
348  /*
349  G4double scale = 1. / G4double(Z);
350  // transform((*prob).begin(),(*prob).end(),(*prob).begin(),bind2nd(multiplies<G4double>(),scale));
351 
352  for (size_t i=0; i<occupancy.size(); i++)
353  {
354  (*prob)[i] *= scale;
355  }
356  */
357  }
358  }
359 }
360 
361 
363 {
364  if (Z < zMin || Z > zMax) {
365 
366  G4Exception("G4ShellData::SelectrandomShell","de0001",FatalErrorInArgument, "Z outside boundaries");
367 
368  }
369  G4int shellIndex = 0;
370  std::vector<G4double> prob = ShellVector(Z);
371  G4double random = G4UniformRand();
372 
373  // std::vector<G4double>::const_iterator pos;
374  // pos = lower_bound(prob.begin(),prob.end(),random);
375 
376  // Binary search the shell with probability less or equal random
377 
378  G4int nShellsLocal = NumberOfShells(Z);
379  G4int upperBound = nShellsLocal;
380 
381  while (shellIndex <= upperBound)
382  {
383  G4int midShell = (shellIndex + upperBound) / 2;
384  if ( random < prob[midShell] )
385  upperBound = midShell - 1;
386  else
387  shellIndex = midShell + 1;
388  }
389  if (shellIndex >= nShellsLocal) shellIndex = nShellsLocal - 1;
390 
391  return shellIndex;
392 }