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