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