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