| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // 23 //
>> 24 // $Id: G4ShellData.cc,v 1.5 2002/05/28 09:20:21 pia Exp $
>> 25 // GEANT4 tag $Name: geant4-05-00 $
27 // 26 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 27 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 28 //
30 // History: 29 // History:
31 // ----------- 30 // -----------
32 // 31 Jul 2001 MGP Created 31 // 31 Jul 2001 MGP Created
33 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war <<
34 // 32 //
35 // ------------------------------------------- 33 // -------------------------------------------------------------------
36 34
37 #include "G4ShellData.hh" 35 #include "G4ShellData.hh"
38 #include "G4DataVector.hh" 36 #include "G4DataVector.hh"
39 #include "G4SystemOfUnits.hh" << 37 #include "g4std/fstream"
40 #include <fstream> << 38 #include "g4std/strstream"
41 #include <sstream> <<
42 #include <numeric> <<
43 #include <algorithm> <<
44 #include <valarray> <<
45 #include <functional> <<
46 #include "Randomize.hh" <<
47 <<
48 //....oooOO0OOooo........oooOO0OOooo........oo <<
49 <<
50 G4ShellData::G4ShellData(G4int minZ, G4int max <<
51 : zMin(minZ), zMax(maxZ), occupancyData(isOc <<
52 {} <<
53 39
54 //....oooOO0OOooo........oooOO0OOooo........oo << 40 // Constructor
55 41
>> 42 G4ShellData::G4ShellData(G4int minZ, G4int maxZ)
>> 43 : zMin(minZ), zMax(maxZ)
>> 44 { }
>> 45
>> 46 // Destructor
56 G4ShellData::~G4ShellData() 47 G4ShellData::~G4ShellData()
57 { 48 {
58 for (auto& pos : idMap) << 49 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::iterator pos;
>> 50
>> 51 for (pos = idMap.begin(); pos != idMap.end(); ++pos)
59 { 52 {
60 std::vector<G4double>* dataSet = pos.sec << 53 G4DataVector* dataSet = (*pos).second;
61 delete dataSet; 54 delete dataSet;
62 } 55 }
63 for (auto& pos2 : bindingMap) << 56 for (pos = bindingMap.begin(); pos != bindingMap.end(); ++pos)
64 { 57 {
65 G4DataVector* dataSet = pos2.second; << 58 G4DataVector* dataSet = (*pos).second;
66 delete dataSet; 59 delete dataSet;
67 } 60 }
68 <<
69 if (occupancyData) <<
70 { <<
71 for (auto& pos3 : occupancyPdfMap) <<
72 { <<
73 std::vector<G4double>* dataSet = pos3.seco <<
74 delete dataSet; <<
75 } <<
76 } <<
77 } 61 }
78 62
79 //....oooOO0OOooo........oooOO0OOooo........oo <<
80 63
81 std::size_t G4ShellData::NumberOfShells(G4int << 64 size_t G4ShellData::NumberOfShells(G4int Z) const
82 { 65 {
83 G4int z = Z - 1; 66 G4int z = Z - 1;
84 G4int n = 0; 67 G4int n = 0;
85 68
86 if (Z>= zMin && Z <= zMax) 69 if (Z>= zMin && Z <= zMax)
87 { 70 {
88 n = nShells[z]; 71 n = nShells[z];
89 } 72 }
90 return n; 73 return n;
91 } 74 }
92 75
93 //....oooOO0OOooo........oooOO0OOooo........oo <<
94 76
95 const std::vector<G4double>& G4ShellData::Shel << 77 const G4DataVector& G4ShellData::ShellIdVector(G4int Z) const
96 { 78 {
97 if (Z < zMin || Z > zMax) { << 79 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::const_iterator pos;
98 G4Exception("G4ShellData::ShellIdVector"," << 80 if (Z < zMin || Z > zMax)
99 } << 81 G4Exception("G4ShellData::ShellIdVector - Z outside boundaries");
100 auto pos = idMap.find(Z); << 82 pos = idMap.find(Z);
101 std::vector<G4double>* dataSet = (*pos).seco << 83 G4DataVector* dataSet = (*pos).second;
102 return *dataSet; 84 return *dataSet;
103 } 85 }
104 86
105 //....oooOO0OOooo........oooOO0OOooo........oo <<
106 <<
107 const std::vector<G4double>& G4ShellData::Shel <<
108 { <<
109 if (Z < zMin || Z > zMax) <<
110 G4Exception("G4ShellData::ShellVector()"," <<
111 auto pos = occupancyPdfMap.find(Z); <<
112 std::vector<G4double>* dataSet = (*pos).seco <<
113 return *dataSet; <<
114 } <<
115 <<
116 //....oooOO0OOooo........oooOO0OOooo........oo <<
117 <<
118 G4int G4ShellData::ShellId(G4int Z, G4int shel 87 G4int G4ShellData::ShellId(G4int Z, G4int shellIndex) const
119 { 88 {
120 G4int n = -1; 89 G4int n = -1;
121 90
122 if (Z >= zMin && Z <= zMax) 91 if (Z >= zMin && Z <= zMax)
123 { 92 {
124 auto pos = idMap.find(Z); << 93 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::const_iterator pos;
>> 94 pos = idMap.find(Z);
125 if (pos!= idMap.end()) 95 if (pos!= idMap.end())
126 { 96 {
127 std::vector<G4double> dataSet = *((*pos).s << 97 G4DataVector dataSet = *((*pos).second);
128 G4int nData = (G4int)dataSet.size(); << 98 G4int nData = dataSet.size();
129 if (shellIndex >= 0 && shellIndex < nData) 99 if (shellIndex >= 0 && shellIndex < nData)
130 { 100 {
131 n = (G4int) dataSet[shellIndex]; 101 n = (G4int) dataSet[shellIndex];
132 } 102 }
133 } 103 }
134 } 104 }
135 return n; 105 return n;
136 } 106 }
137 107
138 //....oooOO0OOooo........oooOO0OOooo........oo <<
139 <<
140 G4double G4ShellData::ShellOccupancyProbabilit <<
141 { <<
142 G4double prob = -1.; <<
143 if (Z >= zMin && Z <= zMax) <<
144 { <<
145 auto pos = idMap.find(Z); <<
146 if (pos!= idMap.end()) <<
147 { <<
148 std::vector<G4double> dataSet = *((*pos).s <<
149 G4int nData = (G4int)dataSet.size(); <<
150 if (shellIndex >= 0 && shellIndex < nData) <<
151 { <<
152 prob = dataSet[shellIndex]; <<
153 } <<
154 } <<
155 } <<
156 return prob; <<
157 } <<
158 <<
159 //....oooOO0OOooo........oooOO0OOooo........oo <<
160 108
161 G4double G4ShellData::BindingEnergy(G4int Z, G 109 G4double G4ShellData::BindingEnergy(G4int Z, G4int shellIndex) const
162 { 110 {
163 G4double value = 0.; 111 G4double value = 0.;
164 112
165 if (Z >= zMin && Z <= zMax) 113 if (Z >= zMin && Z <= zMax)
166 { 114 {
167 auto pos = bindingMap.find(Z); << 115 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::const_iterator pos;
>> 116 pos = bindingMap.find(Z);
168 if (pos!= bindingMap.end()) 117 if (pos!= bindingMap.end())
169 { 118 {
170 G4DataVector dataSet = *((*pos).second); 119 G4DataVector dataSet = *((*pos).second);
171 G4int nData = (G4int)dataSet.size(); << 120 G4int nData = dataSet.size();
172 if (shellIndex >= 0 && shellIndex < nData) 121 if (shellIndex >= 0 && shellIndex < nData)
173 { 122 {
174 value = dataSet[shellIndex]; 123 value = dataSet[shellIndex];
175 } 124 }
176 } 125 }
177 } 126 }
178 return value; 127 return value;
179 } 128 }
180 129
181 //....oooOO0OOooo........oooOO0OOooo........oo <<
182 <<
183 void G4ShellData::PrintData() const 130 void G4ShellData::PrintData() const
184 { 131 {
185 for (G4int Z = zMin; Z <= zMax; Z++) 132 for (G4int Z = zMin; Z <= zMax; Z++)
186 { 133 {
187 G4cout << "---- Shell data for Z = " 134 G4cout << "---- Shell data for Z = "
188 << Z 135 << Z
189 << " ---- " 136 << " ---- "
190 << G4endl; 137 << G4endl;
191 G4int nSh = nShells[Z-1]; 138 G4int nSh = nShells[Z-1];
192 auto posId = idMap.find(Z); << 139 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::const_iterator posId;
193 std::vector<G4double>* ids = (*posId).se << 140 posId = idMap.find(Z);
194 auto posE = bindingMap.find(Z); << 141 G4DataVector* ids = (*posId).second;
>> 142 G4std::map<G4int,G4DataVector*,G4std::less<G4int> >::const_iterator posE;
>> 143 posE = bindingMap.find(Z);
195 G4DataVector* energies = (*posE).second; 144 G4DataVector* energies = (*posE).second;
196 for (G4int i=0; i<nSh; ++i) << 145 for (G4int i=0; i<nSh; i++)
197 { 146 {
198 G4int id = (G4int) (*ids)[i]; 147 G4int id = (G4int) (*ids)[i];
199 G4double e = (*energies)[i] / keV; << 148 G4double e = (*energies)[i] / MeV;
200 G4cout << i << ") "; << 149 G4cout << i <<") Shell id: " << id
201 << 150 << " - Binding energy = "
202 if (occupancyData) << 151 << e << " MeV " << G4endl;
203 { <<
204 G4cout << " Occupancy: "; <<
205 } <<
206 else <<
207 { <<
208 G4cout << " Shell id: "; <<
209 } <<
210 G4cout << id << " - Binding energy = " <<
211 << e << " keV "; <<
212 if (occupancyData) <<
213 { <<
214 auto posOcc = occupancyPdfMap.find(Z); <<
215 std::vector<G4double> probs = <<
216 G4double prob = probs[i]; <<
217 G4cout << "- Probability = " << prob; <<
218 } <<
219 G4cout << G4endl; <<
220 } 152 }
221 G4cout << "----------------------------- 153 G4cout << "-------------------------------------------------"
222 << G4endl; 154 << G4endl;
223 } 155 }
224 } 156 }
225 157
226 //....oooOO0OOooo........oooOO0OOooo........oo <<
227 158
228 void G4ShellData::LoadData(const G4String& fil 159 void G4ShellData::LoadData(const G4String& fileName)
229 { 160 {
230 // Build the complete string identifying the 161 // Build the complete string identifying the file with the data set
231 std::ostringstream ost; << 162
>> 163 char nameChar[100] = {""};
>> 164 G4std::ostrstream ost(nameChar, 100, G4std::ios::out);
>> 165
232 ost << fileName << ".dat"; 166 ost << fileName << ".dat";
233 G4String name(ost.str()); <<
234 167
235 const char* path = G4FindDataDir("G4LEDATA") << 168 G4String name(nameChar);
>> 169
>> 170 char* path = getenv("G4LEDATA");
236 if (!path) 171 if (!path)
237 { 172 {
238 G4String excep("G4ShellData::LoadData()" << 173 G4String excep("G4EMDataSet - G4LEDATA environment variable not set");
239 G4Exception(excep,"em0006",FatalExceptio << 174 G4Exception(excep);
240 return; <<
241 } 175 }
242 176
243 G4String pathString(path); 177 G4String pathString(path);
244 G4String dirFile = pathString + name; 178 G4String dirFile = pathString + name;
245 std::ifstream file(dirFile); << 179 G4std::ifstream file(dirFile);
246 std::filebuf* lsdp = file.rdbuf(); << 180 G4std::filebuf* lsdp = file.rdbuf();
247 181
248 if (! (lsdp->is_open()) ) 182 if (! (lsdp->is_open()) )
249 { 183 {
250 << 184 G4String s1("G4ShellData - data file: ");
251 G4String excep = "G4ShellData::LoadData( << 185 G4String s2(" not found");
252 G4String msg = "data file: " + dirFile + << 186 G4String excep = s1 + dirFile + s2;
253 G4Exception(excep, "em0003",FatalExcepti << 187 G4Exception(excep);
254 return; <<
255 } 188 }
256 189
257 G4double a = 0; 190 G4double a = 0;
258 G4int k = 1; 191 G4int k = 1;
259 G4int sLocal = 0; << 192 G4int s = 0;
260 193
261 G4int Z = 1; 194 G4int Z = 1;
262 G4DataVector* energies = new G4DataVector; 195 G4DataVector* energies = new G4DataVector;
263 std::vector<G4double>* ids = new std::vector << 196 G4DataVector* ids = new G4DataVector;
264 197
265 do { 198 do {
266 file >> a; 199 file >> a;
267 G4int nColumns = 2; 200 G4int nColumns = 2;
268 if (a == -1) 201 if (a == -1)
269 { 202 {
270 if (sLocal == 0) << 203 if (s == 0)
271 { 204 {
272 // End of a shell data set 205 // End of a shell data set
273 idMap[Z] = ids; 206 idMap[Z] = ids;
274 bindingMap[Z] = energies; 207 bindingMap[Z] = energies;
275 G4int n = (G4int)ids->size(); << 208 G4int n = ids->size();
276 nShells.push_back(n); 209 nShells.push_back(n);
277 // Start of new shell data set 210 // Start of new shell data set
278 << 211 ids = new G4DataVector;
279 ids = new std::vector<G4double>; <<
280 energies = new G4DataVector; 212 energies = new G4DataVector;
281 Z++; 213 Z++;
282 } 214 }
283 sLocal++; << 215 s++;
284 if (sLocal == nColumns) << 216 if (s == nColumns)
285 { 217 {
286 sLocal = 0; << 218 s = 0;
287 } 219 }
288 } 220 }
289 << 221 else if (a == -2)
290 // moved out of the do-while since might g << 222 {
291 // else if (a == -2) << 223 // End of file; delete the empty vectors created when encountering the last -1 -1 row
292 // { << 224 delete energies;
293 // End of file; delete the empty vectors c << 225 delete ids;
294 // delete energies; << 226 //nComponents = components.size();
295 // delete ids; << 227 }
296 //nComponents = components.size(); <<
297 // } <<
298 else 228 else
299 { 229 {
300 // 1st column is shell id 230 // 1st column is shell id
301 if(k%nColumns != 0) 231 if(k%nColumns != 0)
302 { 232 {
303 ids->push_back(a); 233 ids->push_back(a);
304 k++; 234 k++;
305 } 235 }
306 else if (k%nColumns == 0) 236 else if (k%nColumns == 0)
307 { 237 {
308 // 2nd column is binding energy 238 // 2nd column is binding energy
309 G4double e = a * MeV; 239 G4double e = a * MeV;
310 energies->push_back(e); 240 energies->push_back(e);
311 k = 1; 241 k = 1;
312 } 242 }
313 } 243 }
314 } while (a != -2); // end of file 244 } while (a != -2); // end of file
315 file.close(); 245 file.close();
316 delete energies; <<
317 delete ids; <<
318 <<
319 // For Doppler broadening: the data set cont <<
320 // Build additional map with probability for <<
321 <<
322 if (occupancyData) <<
323 { <<
324 // Build cumulative from raw shell occup <<
325 for (G4int ZLocal=zMin; ZLocal <= zMax; <<
326 { <<
327 std::vector<G4double> occupancy = ShellIdV <<
328 <<
329 std::vector<G4double>* prob = new std::vec <<
330 G4double scale = 1. / G4double(ZLocal); <<
331 <<
332 prob->push_back(occupancy[0] * scale); <<
333 for (std::size_t i=1; i<occupancy.size(); <<
334 { <<
335 prob->push_back(occupancy[i]*scale + ( <<
336 } <<
337 occupancyPdfMap[ZLocal] = prob; <<
338 } <<
339 } <<
340 } 246 }
341 247
342 //....oooOO0OOooo........oooOO0OOooo........oo <<
343 <<
344 G4int G4ShellData::SelectRandomShell(G4int Z) <<
345 { <<
346 if (Z < zMin || Z > zMax) <<
347 G4Exception("G4ShellData::SelectrandomShel <<
348 <<
349 G4int shellIndex = 0; <<
350 std::vector<G4double> prob = ShellVector(Z); <<
351 G4double random = G4UniformRand(); <<
352 <<
353 // Binary search the shell with probability <<
354 G4int nShellsLocal = (G4int)NumberOfShells(Z <<
355 G4int upperBound = nShellsLocal; <<
356 <<
357 while (shellIndex <= upperBound) <<
358 { <<
359 G4int midShell = (shellIndex + upperBoun <<
360 if ( random < prob[midShell] ) <<
361 upperBound = midShell - 1; <<
362 else <<
363 shellIndex = midShell + 1; <<
364 } <<
365 if (shellIndex >= nShellsLocal) shellIndex = <<
366 <<
367 return shellIndex; <<
368 } <<
369 248