| 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.12 2010-12-27 17:44:50 vnivanch Exp $
>> 28 // GEANT4 tag $Name: not supported by cvs2svn $
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 35 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings
34 // 36 //
35 // ------------------------------------------- 37 // -------------------------------------------------------------------
36 38
37 #include "G4ShellData.hh" 39 #include "G4ShellData.hh"
38 #include "G4DataVector.hh" 40 #include "G4DataVector.hh"
39 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh"
40 #include <fstream> 42 #include <fstream>
41 #include <sstream> 43 #include <sstream>
42 #include <numeric> 44 #include <numeric>
43 #include <algorithm> 45 #include <algorithm>
44 #include <valarray> 46 #include <valarray>
45 #include <functional> 47 #include <functional>
46 #include "Randomize.hh" 48 #include "Randomize.hh"
47 49
48 //....oooOO0OOooo........oooOO0OOooo........oo << 50 // Constructor
49 51
50 G4ShellData::G4ShellData(G4int minZ, G4int max 52 G4ShellData::G4ShellData(G4int minZ, G4int maxZ, G4bool isOccupancy)
51 : zMin(minZ), zMax(maxZ), occupancyData(isOc 53 : zMin(minZ), zMax(maxZ), occupancyData(isOccupancy)
52 {} << 54 { }
53 <<
54 //....oooOO0OOooo........oooOO0OOooo........oo <<
55 55
>> 56 // Destructor
56 G4ShellData::~G4ShellData() 57 G4ShellData::~G4ShellData()
57 { 58 {
58 for (auto& pos : idMap) << 59 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos;
>> 60 for (pos = idMap.begin(); pos != idMap.end(); ++pos)
59 { 61 {
60 std::vector<G4double>* dataSet = pos.sec << 62 std::vector<G4double>* dataSet = (*pos).second;
61 delete dataSet; 63 delete dataSet;
62 } 64 }
63 for (auto& pos2 : bindingMap) << 65
>> 66 std::map<G4int,G4DataVector*,std::less<G4int> >::iterator pos2;
>> 67 for (pos2 = bindingMap.begin(); pos2 != bindingMap.end(); ++pos2)
64 { 68 {
65 G4DataVector* dataSet = pos2.second; << 69 G4DataVector* dataSet = (*pos2).second;
66 delete dataSet; 70 delete dataSet;
67 } 71 }
68 72
69 if (occupancyData) 73 if (occupancyData)
70 { 74 {
71 for (auto& pos3 : occupancyPdfMap) << 75 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::iterator pos3;
>> 76 for (pos3 = occupancyPdfMap.begin(); pos3 != occupancyPdfMap.end(); ++pos3)
72 { 77 {
73 std::vector<G4double>* dataSet = pos3.seco << 78 std::vector<G4double>* dataSet = (*pos3).second;
74 delete dataSet; 79 delete dataSet;
75 } 80 }
76 } 81 }
77 } 82 }
78 83
79 //....oooOO0OOooo........oooOO0OOooo........oo <<
80 84
81 std::size_t G4ShellData::NumberOfShells(G4int << 85 size_t G4ShellData::NumberOfShells(G4int Z) const
82 { 86 {
83 G4int z = Z - 1; 87 G4int z = Z - 1;
84 G4int n = 0; 88 G4int n = 0;
85 89
86 if (Z>= zMin && Z <= zMax) 90 if (Z>= zMin && Z <= zMax)
87 { 91 {
88 n = nShells[z]; 92 n = nShells[z];
89 } 93 }
90 return n; 94 return n;
91 } 95 }
92 96
93 //....oooOO0OOooo........oooOO0OOooo........oo <<
94 97
95 const std::vector<G4double>& G4ShellData::Shel 98 const std::vector<G4double>& G4ShellData::ShellIdVector(G4int Z) const
96 { 99 {
97 if (Z < zMin || Z > zMax) { << 100 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
98 G4Exception("G4ShellData::ShellIdVector"," << 101 if (Z < zMin || Z > zMax) G4Exception("G4ShellData::ShellIdVector - Z outside boundaries");
99 } << 102 pos = idMap.find(Z);
100 auto pos = idMap.find(Z); <<
101 std::vector<G4double>* dataSet = (*pos).seco 103 std::vector<G4double>* dataSet = (*pos).second;
102 return *dataSet; 104 return *dataSet;
103 } 105 }
104 106
105 //....oooOO0OOooo........oooOO0OOooo........oo <<
106 107
107 const std::vector<G4double>& G4ShellData::Shel 108 const std::vector<G4double>& G4ShellData::ShellVector(G4int Z) const
108 { 109 {
109 if (Z < zMin || Z > zMax) << 110 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
110 G4Exception("G4ShellData::ShellVector()"," << 111 if (Z < zMin || Z > zMax) G4Exception("G4ShellData::ShellVector - Z outside boundaries");
111 auto pos = occupancyPdfMap.find(Z); << 112 pos = occupancyPdfMap.find(Z);
112 std::vector<G4double>* dataSet = (*pos).seco 113 std::vector<G4double>* dataSet = (*pos).second;
113 return *dataSet; 114 return *dataSet;
114 } 115 }
115 116
116 //....oooOO0OOooo........oooOO0OOooo........oo <<
117 117
118 G4int G4ShellData::ShellId(G4int Z, G4int shel 118 G4int G4ShellData::ShellId(G4int Z, G4int shellIndex) const
119 { 119 {
120 G4int n = -1; 120 G4int n = -1;
121 121
122 if (Z >= zMin && Z <= zMax) 122 if (Z >= zMin && Z <= zMax)
123 { 123 {
124 auto pos = idMap.find(Z); << 124 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
>> 125 pos = idMap.find(Z);
125 if (pos!= idMap.end()) 126 if (pos!= idMap.end())
126 { 127 {
127 std::vector<G4double> dataSet = *((*pos).s 128 std::vector<G4double> dataSet = *((*pos).second);
128 G4int nData = (G4int)dataSet.size(); << 129 G4int nData = dataSet.size();
129 if (shellIndex >= 0 && shellIndex < nData) 130 if (shellIndex >= 0 && shellIndex < nData)
130 { 131 {
131 n = (G4int) dataSet[shellIndex]; 132 n = (G4int) dataSet[shellIndex];
132 } 133 }
133 } 134 }
134 } 135 }
135 return n; 136 return n;
136 } 137 }
137 138
138 //....oooOO0OOooo........oooOO0OOooo........oo <<
139 139
140 G4double G4ShellData::ShellOccupancyProbabilit 140 G4double G4ShellData::ShellOccupancyProbability(G4int Z, G4int shellIndex) const
141 { 141 {
142 G4double prob = -1.; 142 G4double prob = -1.;
>> 143
143 if (Z >= zMin && Z <= zMax) 144 if (Z >= zMin && Z <= zMax)
144 { 145 {
145 auto pos = idMap.find(Z); << 146 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator pos;
>> 147 pos = idMap.find(Z);
146 if (pos!= idMap.end()) 148 if (pos!= idMap.end())
147 { 149 {
148 std::vector<G4double> dataSet = *((*pos).s 150 std::vector<G4double> dataSet = *((*pos).second);
149 G4int nData = (G4int)dataSet.size(); << 151 G4int nData = dataSet.size();
150 if (shellIndex >= 0 && shellIndex < nData) 152 if (shellIndex >= 0 && shellIndex < nData)
151 { 153 {
152 prob = dataSet[shellIndex]; 154 prob = dataSet[shellIndex];
153 } 155 }
154 } 156 }
155 } 157 }
156 return prob; 158 return prob;
157 } 159 }
158 160
159 //....oooOO0OOooo........oooOO0OOooo........oo << 161
160 162
161 G4double G4ShellData::BindingEnergy(G4int Z, G 163 G4double G4ShellData::BindingEnergy(G4int Z, G4int shellIndex) const
162 { 164 {
163 G4double value = 0.; 165 G4double value = 0.;
164 166
165 if (Z >= zMin && Z <= zMax) 167 if (Z >= zMin && Z <= zMax)
166 { 168 {
167 auto pos = bindingMap.find(Z); << 169 std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator pos;
>> 170 pos = bindingMap.find(Z);
168 if (pos!= bindingMap.end()) 171 if (pos!= bindingMap.end())
169 { 172 {
170 G4DataVector dataSet = *((*pos).second); 173 G4DataVector dataSet = *((*pos).second);
171 G4int nData = (G4int)dataSet.size(); << 174 G4int nData = dataSet.size();
172 if (shellIndex >= 0 && shellIndex < nData) 175 if (shellIndex >= 0 && shellIndex < nData)
173 { 176 {
174 value = dataSet[shellIndex]; 177 value = dataSet[shellIndex];
175 } 178 }
176 } 179 }
177 } 180 }
178 return value; 181 return value;
179 } 182 }
180 183
181 //....oooOO0OOooo........oooOO0OOooo........oo <<
182 <<
183 void G4ShellData::PrintData() const 184 void G4ShellData::PrintData() const
184 { 185 {
185 for (G4int Z = zMin; Z <= zMax; Z++) 186 for (G4int Z = zMin; Z <= zMax; Z++)
186 { 187 {
187 G4cout << "---- Shell data for Z = " 188 G4cout << "---- Shell data for Z = "
188 << Z 189 << Z
189 << " ---- " 190 << " ---- "
190 << G4endl; 191 << G4endl;
191 G4int nSh = nShells[Z-1]; 192 G4int nSh = nShells[Z-1];
192 auto posId = idMap.find(Z); << 193 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posId;
>> 194 posId = idMap.find(Z);
193 std::vector<G4double>* ids = (*posId).se 195 std::vector<G4double>* ids = (*posId).second;
194 auto posE = bindingMap.find(Z); << 196 std::map<G4int,G4DataVector*,std::less<G4int> >::const_iterator posE;
>> 197 posE = bindingMap.find(Z);
195 G4DataVector* energies = (*posE).second; 198 G4DataVector* energies = (*posE).second;
196 for (G4int i=0; i<nSh; ++i) << 199 for (G4int i=0; i<nSh; i++)
197 { 200 {
198 G4int id = (G4int) (*ids)[i]; 201 G4int id = (G4int) (*ids)[i];
199 G4double e = (*energies)[i] / keV; 202 G4double e = (*energies)[i] / keV;
200 G4cout << i << ") "; 203 G4cout << i << ") ";
201 204
202 if (occupancyData) 205 if (occupancyData)
203 { 206 {
204 G4cout << " Occupancy: "; 207 G4cout << " Occupancy: ";
205 } 208 }
206 else 209 else
207 { 210 {
208 G4cout << " Shell id: "; 211 G4cout << " Shell id: ";
209 } 212 }
210 G4cout << id << " - Binding energy = " 213 G4cout << id << " - Binding energy = "
211 << e << " keV "; 214 << e << " keV ";
212 if (occupancyData) 215 if (occupancyData)
213 { 216 {
214 auto posOcc = occupancyPdfMap.find(Z); << 217 std::map<G4int,std::vector<G4double>*,std::less<G4int> >::const_iterator posOcc;
>> 218 posOcc = occupancyPdfMap.find(Z);
215 std::vector<G4double> probs = 219 std::vector<G4double> probs = *((*posOcc).second);
216 G4double prob = probs[i]; 220 G4double prob = probs[i];
217 G4cout << "- Probability = " << prob; 221 G4cout << "- Probability = " << prob;
218 } 222 }
219 G4cout << G4endl; 223 G4cout << G4endl;
220 } 224 }
221 G4cout << "----------------------------- 225 G4cout << "-------------------------------------------------"
222 << G4endl; 226 << G4endl;
223 } 227 }
224 } 228 }
225 229
226 //....oooOO0OOooo........oooOO0OOooo........oo <<
227 230
228 void G4ShellData::LoadData(const G4String& fil 231 void G4ShellData::LoadData(const G4String& fileName)
229 { 232 {
230 // Build the complete string identifying the 233 // Build the complete string identifying the file with the data set
>> 234
231 std::ostringstream ost; 235 std::ostringstream ost;
>> 236
232 ost << fileName << ".dat"; 237 ost << fileName << ".dat";
>> 238
233 G4String name(ost.str()); 239 G4String name(ost.str());
234 240
235 const char* path = G4FindDataDir("G4LEDATA") << 241 char* path = getenv("G4LEDATA");
236 if (!path) 242 if (!path)
237 { 243 {
238 G4String excep("G4ShellData::LoadData()" << 244 G4String excep("G4EMDataSet - G4LEDATA environment variable not set");
239 G4Exception(excep,"em0006",FatalExceptio << 245 G4Exception(excep);
240 return; 246 return;
241 } 247 }
242 248
243 G4String pathString(path); 249 G4String pathString(path);
244 G4String dirFile = pathString + name; 250 G4String dirFile = pathString + name;
245 std::ifstream file(dirFile); 251 std::ifstream file(dirFile);
246 std::filebuf* lsdp = file.rdbuf(); 252 std::filebuf* lsdp = file.rdbuf();
247 253
248 if (! (lsdp->is_open()) ) 254 if (! (lsdp->is_open()) )
249 { 255 {
250 << 256 G4String s1("G4ShellData - data file: ");
251 G4String excep = "G4ShellData::LoadData( << 257 G4String s2(" not found");
252 G4String msg = "data file: " + dirFile + << 258 G4String excep = s1 + dirFile + s2;
253 G4Exception(excep, "em0003",FatalExcepti << 259 G4Exception(excep);
254 return; 260 return;
255 } 261 }
256 262
257 G4double a = 0; 263 G4double a = 0;
258 G4int k = 1; 264 G4int k = 1;
259 G4int sLocal = 0; << 265 G4int s = 0;
260 266
261 G4int Z = 1; 267 G4int Z = 1;
262 G4DataVector* energies = new G4DataVector; 268 G4DataVector* energies = new G4DataVector;
263 std::vector<G4double>* ids = new std::vector 269 std::vector<G4double>* ids = new std::vector<G4double>;
264 270
265 do { 271 do {
266 file >> a; 272 file >> a;
267 G4int nColumns = 2; 273 G4int nColumns = 2;
268 if (a == -1) 274 if (a == -1)
269 { 275 {
270 if (sLocal == 0) << 276 if (s == 0)
271 { 277 {
272 // End of a shell data set 278 // End of a shell data set
273 idMap[Z] = ids; 279 idMap[Z] = ids;
274 bindingMap[Z] = energies; 280 bindingMap[Z] = energies;
275 G4int n = (G4int)ids->size(); << 281 G4int n = ids->size();
276 nShells.push_back(n); 282 nShells.push_back(n);
277 // Start of new shell data set 283 // Start of new shell data set
278 <<
279 ids = new std::vector<G4double>; 284 ids = new std::vector<G4double>;
280 energies = new G4DataVector; 285 energies = new G4DataVector;
281 Z++; 286 Z++;
282 } 287 }
283 sLocal++; << 288 s++;
284 if (sLocal == nColumns) << 289 if (s == nColumns)
285 { 290 {
286 sLocal = 0; << 291 s = 0;
287 } 292 }
288 } 293 }
289 << 294 else if (a == -2)
290 // moved out of the do-while since might g << 295 {
291 // else if (a == -2) << 296 // End of file; delete the empty vectors created when encountering the last -1 -1 row
292 // { << 297 delete energies;
293 // End of file; delete the empty vectors c << 298 delete ids;
294 // delete energies; << 299 //nComponents = components.size();
295 // delete ids; << 300 }
296 //nComponents = components.size(); <<
297 // } <<
298 else 301 else
299 { 302 {
300 // 1st column is shell id 303 // 1st column is shell id
301 if(k%nColumns != 0) 304 if(k%nColumns != 0)
302 { 305 {
303 ids->push_back(a); 306 ids->push_back(a);
304 k++; 307 k++;
305 } 308 }
306 else if (k%nColumns == 0) 309 else if (k%nColumns == 0)
307 { 310 {
308 // 2nd column is binding energy 311 // 2nd column is binding energy
309 G4double e = a * MeV; 312 G4double e = a * MeV;
310 energies->push_back(e); 313 energies->push_back(e);
311 k = 1; 314 k = 1;
312 } 315 }
313 } 316 }
314 } while (a != -2); // end of file 317 } while (a != -2); // end of file
315 file.close(); 318 file.close();
316 delete energies; <<
317 delete ids; <<
318 319
319 // For Doppler broadening: the data set cont 320 // For Doppler broadening: the data set contains shell occupancy and binding energy for each shell
320 // Build additional map with probability for 321 // Build additional map with probability for each shell based on its occupancy
321 322
322 if (occupancyData) 323 if (occupancyData)
323 { 324 {
324 // Build cumulative from raw shell occup 325 // Build cumulative from raw shell occupancy
325 for (G4int ZLocal=zMin; ZLocal <= zMax; << 326
>> 327 for (G4int Z=zMin; Z <= zMax; Z++)
326 { 328 {
327 std::vector<G4double> occupancy = ShellIdV << 329 std::vector<G4double> occupancy = ShellIdVector(Z);
328 330
329 std::vector<G4double>* prob = new std::vec 331 std::vector<G4double>* prob = new std::vector<G4double>;
330 G4double scale = 1. / G4double(ZLocal); << 332 G4double scale = 1. / G4double(Z);
331 333
332 prob->push_back(occupancy[0] * scale); 334 prob->push_back(occupancy[0] * scale);
333 for (std::size_t i=1; i<occupancy.size(); << 335 for (size_t i=1; i<occupancy.size(); i++)
334 { 336 {
335 prob->push_back(occupancy[i]*scale + ( 337 prob->push_back(occupancy[i]*scale + (*prob)[i-1]);
336 } 338 }
337 occupancyPdfMap[ZLocal] = prob; << 339 occupancyPdfMap[Z] = prob;
>> 340
>> 341 /*
>> 342 G4double scale = 1. / G4double(Z);
>> 343 // transform((*prob).begin(),(*prob).end(),(*prob).begin(),bind2nd(multiplies<G4double>(),scale));
>> 344
>> 345 for (size_t i=0; i<occupancy.size(); i++)
>> 346 {
>> 347 (*prob)[i] *= scale;
>> 348 }
>> 349 */
338 } 350 }
339 } 351 }
340 } 352 }
341 353
342 //....oooOO0OOooo........oooOO0OOooo........oo <<
343 354
344 G4int G4ShellData::SelectRandomShell(G4int Z) 355 G4int G4ShellData::SelectRandomShell(G4int Z) const
345 { 356 {
346 if (Z < zMin || Z > zMax) << 357 if (Z < zMin || Z > zMax) G4Exception("G4ShellData::RandomSelect - Z outside boundaries");
347 G4Exception("G4ShellData::SelectrandomShel <<
348 358
349 G4int shellIndex = 0; 359 G4int shellIndex = 0;
350 std::vector<G4double> prob = ShellVector(Z); 360 std::vector<G4double> prob = ShellVector(Z);
351 G4double random = G4UniformRand(); 361 G4double random = G4UniformRand();
352 362
>> 363 // std::vector<G4double>::const_iterator pos;
>> 364 // pos = lower_bound(prob.begin(),prob.end(),random);
>> 365
353 // Binary search the shell with probability 366 // Binary search the shell with probability less or equal random
354 G4int nShellsLocal = (G4int)NumberOfShells(Z << 367
355 G4int upperBound = nShellsLocal; << 368 G4int nShells = NumberOfShells(Z);
>> 369 G4int upperBound = nShells;
356 370
357 while (shellIndex <= upperBound) 371 while (shellIndex <= upperBound)
358 { 372 {
359 G4int midShell = (shellIndex + upperBoun 373 G4int midShell = (shellIndex + upperBound) / 2;
360 if ( random < prob[midShell] ) 374 if ( random < prob[midShell] )
361 upperBound = midShell - 1; 375 upperBound = midShell - 1;
362 else 376 else
363 shellIndex = midShell + 1; 377 shellIndex = midShell + 1;
364 } 378 }
365 if (shellIndex >= nShellsLocal) shellIndex = << 379 if (shellIndex >= nShells) shellIndex = nShells - 1;
366 380
367 return shellIndex; 381 return shellIndex;
368 } 382 }
369 383