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