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