| 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 // History: 26 // History:
27 // ----------- 27 // -----------
28 // 31 Jul 2001 MGP Created 28 // 31 Jul 2001 MGP Created
29 // 29 //
30 // 15 Jul 2009 Nicolas A. Karakatsanis 30 // 15 Jul 2009 Nicolas A. Karakatsanis
31 // 31 //
32 // - New Constructor 32 // - New Constructor was added when logarithmic data are loaded as well
33 // to enhance CPU 33 // to enhance CPU performance
34 // 34 //
35 // - Destructor was 35 // - Destructor was modified accordingly
36 // 36 //
37 // - LoadNonLogData 37 // - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW
38 // dataset. It is 38 // dataset. It is essentially performing the data loading operations as in the past.
39 // 39 //
40 // - LoadData method 40 // - LoadData method was revised in order to calculate the logarithmic values of the data
41 // It retrieves th 41 // It retrieves the data values from the G4EMLOW data files but, then, calculates the
42 // respective log 42 // respective log values and loads them to seperate data structures.
43 // 43 //
44 // - SetLogEnergiesD 44 // - SetLogEnergiesData method was cretaed to set logarithmic values to G4 data vectors.
45 // The EM data set 45 // The EM data sets, initialized this way, contain both non-log and log values.
46 // These initializ 46 // These initialized data sets can enhance the computing performance of data interpolation
47 // operations 47 // operations
48 // 48 //
49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war 49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings and cleanup logic
50 // 50 //
51 // ------------------------------------------- 51 // -------------------------------------------------------------------
52 52
53 #include "G4EMDataSet.hh" 53 #include "G4EMDataSet.hh"
54 #include "G4VDataSetAlgorithm.hh" 54 #include "G4VDataSetAlgorithm.hh"
55 #include <fstream> 55 #include <fstream>
56 #include <sstream> 56 #include <sstream>
57 #include "G4Integrator.hh" 57 #include "G4Integrator.hh"
58 #include "Randomize.hh" 58 #include "Randomize.hh"
59 #include "G4LinInterpolation.hh" 59 #include "G4LinInterpolation.hh"
60 60
>> 61
61 G4EMDataSet::G4EMDataSet(G4int Z, 62 G4EMDataSet::G4EMDataSet(G4int Z,
62 G4VDataSetAlgorithm* algo, 63 G4VDataSetAlgorithm* algo,
63 G4double xUnit, 64 G4double xUnit,
64 G4double yUnit, 65 G4double yUnit,
65 G4bool random): << 66 G4bool random):
66 energies(nullptr), << 67 z(Z),
67 data(nullptr), << 68 energies(0),
68 log_energies(nullptr), << 69 data(0),
69 log_data(nullptr), << 70 log_energies(0),
>> 71 log_data(0),
70 algorithm(algo), 72 algorithm(algo),
71 pdf(nullptr), <<
72 unitEnergies(xUnit), 73 unitEnergies(xUnit),
73 unitData(yUnit), << 74 unitData(yUnit),
74 z(Z), << 75 pdf(0),
75 randomSet(random) 76 randomSet(random)
76 { 77 {
77 if (algorithm == nullptr) { << 78 if (algorithm == 0) {
78 G4Exception("G4EMDataSet::G4EMDataSet", 79 G4Exception("G4EMDataSet::G4EMDataSet",
79 "em1012",FatalException,"interpolation 80 "em1012",FatalException,"interpolation == 0");
80 } else if (randomSet) { BuildPdf(); } 81 } else if (randomSet) { BuildPdf(); }
81 } 82 }
82 83
83 //....oooOO0OOooo........oooOO0OOooo........oo <<
84 <<
85 G4EMDataSet::G4EMDataSet(G4int argZ, 84 G4EMDataSet::G4EMDataSet(G4int argZ,
86 G4DataVector* dataX, 85 G4DataVector* dataX,
87 G4DataVector* dataY, 86 G4DataVector* dataY,
88 G4VDataSetAlgorithm* algo, 87 G4VDataSetAlgorithm* algo,
89 G4double xUnit, 88 G4double xUnit,
90 G4double yUnit, 89 G4double yUnit,
91 G4bool random): 90 G4bool random):
>> 91 z(argZ),
92 energies(dataX), 92 energies(dataX),
93 data(dataY), 93 data(dataY),
94 log_energies(nullptr), << 94 log_energies(0),
95 log_data(nullptr), << 95 log_data(0),
96 algorithm(algo), 96 algorithm(algo),
97 pdf(nullptr), <<
98 unitEnergies(xUnit), 97 unitEnergies(xUnit),
99 unitData(yUnit), << 98 unitData(yUnit),
100 z(argZ), << 99 pdf(0),
101 randomSet(random) 100 randomSet(random)
102 { 101 {
103 if (!algorithm || !energies || !data) { 102 if (!algorithm || !energies || !data) {
104 G4Exception("G4EMDataSet::G4EMDataSet", 103 G4Exception("G4EMDataSet::G4EMDataSet",
105 "em1012",FatalException,"interpolation 104 "em1012",FatalException,"interpolation == 0");
106 } else { 105 } else {
107 if (energies->size() != data->size()) { 106 if (energies->size() != data->size()) {
108 G4Exception("G4EMDataSet::G4EMDataSet", 107 G4Exception("G4EMDataSet::G4EMDataSet",
109 "em1012",FatalException,"different size 108 "em1012",FatalException,"different size for energies and data");
110 } else if (randomSet) { 109 } else if (randomSet) {
111 BuildPdf(); 110 BuildPdf();
112 } 111 }
113 } 112 }
114 } 113 }
115 114
116 //....oooOO0OOooo........oooOO0OOooo........oo <<
117 <<
118 G4EMDataSet::G4EMDataSet(G4int argZ, 115 G4EMDataSet::G4EMDataSet(G4int argZ,
119 G4DataVector* dataX, 116 G4DataVector* dataX,
120 G4DataVector* dataY, 117 G4DataVector* dataY,
121 G4DataVector* dataLog 118 G4DataVector* dataLogX,
122 G4DataVector* dataLog 119 G4DataVector* dataLogY,
123 G4VDataSetAlgorithm* algo, 120 G4VDataSetAlgorithm* algo,
124 G4double xUnit, 121 G4double xUnit,
125 G4double yUnit, 122 G4double yUnit,
126 G4bool random): 123 G4bool random):
>> 124 z(argZ),
127 energies(dataX), 125 energies(dataX),
128 data(dataY), 126 data(dataY),
129 log_energies(dataLogX), 127 log_energies(dataLogX),
130 log_data(dataLogY), 128 log_data(dataLogY),
131 algorithm(algo), 129 algorithm(algo),
132 pdf(nullptr), <<
133 unitEnergies(xUnit), 130 unitEnergies(xUnit),
134 unitData(yUnit), 131 unitData(yUnit),
135 z(argZ), << 132 pdf(0),
136 randomSet(random) 133 randomSet(random)
137 { 134 {
138 if (!algorithm || !energies || !data || !log 135 if (!algorithm || !energies || !data || !log_energies || !log_data) {
139 G4Exception("G4EMDataSet::G4EMDataSet", 136 G4Exception("G4EMDataSet::G4EMDataSet",
140 "em1012",FatalException,"interpolation 137 "em1012",FatalException,"interpolation == 0");
141 } else { 138 } else {
142 if ((energies->size() != data->size()) || 139 if ((energies->size() != data->size()) ||
143 (energies->size() != log_energies->size()) | 140 (energies->size() != log_energies->size()) ||
144 (energies->size() != log_data->size())) { 141 (energies->size() != log_data->size())) {
145 G4Exception("G4EMDataSet::G4EMDataSet", 142 G4Exception("G4EMDataSet::G4EMDataSet",
146 "em1012",FatalException,"different size 143 "em1012",FatalException,"different size for energies and data");
147 } else if (randomSet) { 144 } else if (randomSet) {
148 BuildPdf(); 145 BuildPdf();
149 } 146 }
150 } 147 }
151 } 148 }
152 149
153 //....oooOO0OOooo........oooOO0OOooo........oo <<
154 <<
155 G4EMDataSet::~G4EMDataSet() 150 G4EMDataSet::~G4EMDataSet()
156 { 151 {
157 delete algorithm; 152 delete algorithm;
158 delete energies; 153 delete energies;
159 delete data; 154 delete data;
160 delete pdf; 155 delete pdf;
161 delete log_energies; 156 delete log_energies;
162 delete log_data; 157 delete log_data;
163 } 158 }
164 159
165 //....oooOO0OOooo........oooOO0OOooo........oo <<
166 <<
167 G4double G4EMDataSet::FindValue(G4double energ 160 G4double G4EMDataSet::FindValue(G4double energy, G4int /* componentId */) const
168 { 161 {
169 if (energy <= (*energies)[0]) { 162 if (energy <= (*energies)[0]) {
170 return (*data)[0]; 163 return (*data)[0];
171 } 164 }
172 std::size_t i = energies->size()-1; << 165 size_t i = energies->size()-1;
173 if (energy >= (*energies)[i]) { return (*dat 166 if (energy >= (*energies)[i]) { return (*data)[i]; }
174 167
175 //Nicolas A. Karakatsanis: Check if the loga 168 //Nicolas A. Karakatsanis: Check if the logarithmic data have been loaded to decide
176 // which Interpolati 169 // which Interpolation-Calculation method will be applied
177 if (log_energies != nullptr) << 170 if (log_energies != 0)
178 { 171 {
179 return algorithm->Calculate(energy, (G4in << 172 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data, *log_energies, *log_data);
180 } 173 }
181 174
182 return algorithm->Calculate(energy, (G4int)F << 175 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
183 } 176 }
184 177
185 //....oooOO0OOooo........oooOO0OOooo........oo <<
186 <<
187 void G4EMDataSet::PrintData(void) const 178 void G4EMDataSet::PrintData(void) const
188 { 179 {
189 std::size_t size = energies->size(); << 180 size_t size = energies->size();
190 for (std::size_t i(0); i<size; ++i) << 181 for (size_t i(0); i<size; i++)
191 { 182 {
192 G4cout << "Point: " << ((*energies)[i]/u 183 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
193 << " - Data value: " << ((*data)[i]/uni 184 << " - Data value: " << ((*data)[i]/unitData);
194 if (pdf != 0) G4cout << " - PDF : " << ( 185 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
195 G4cout << G4endl; 186 G4cout << G4endl;
196 } 187 }
197 } 188 }
198 189
199 //....oooOO0OOooo........oooOO0OOooo........oo <<
200 <<
201 void G4EMDataSet::SetEnergiesData(G4DataVector 190 void G4EMDataSet::SetEnergiesData(G4DataVector* dataX,
202 G4DataVector* dataY, 191 G4DataVector* dataY,
203 G4int /* componentId */) 192 G4int /* componentId */)
204 { 193 {
205 if(!dataX || !dataY) { 194 if(!dataX || !dataY) {
206 G4Exception("G4EMDataSet::SetEnergiesData" 195 G4Exception("G4EMDataSet::SetEnergiesData",
207 "em1012",FatalException,"new interpola 196 "em1012",FatalException,"new interpolation == 0");
208 } else { 197 } else {
209 if (dataX->size() != dataY->size()) { 198 if (dataX->size() != dataY->size()) {
210 G4Exception("G4EMDataSet::SetEnergiesDat 199 G4Exception("G4EMDataSet::SetEnergiesData",
211 "em1012",FatalException,"different size 200 "em1012",FatalException,"different size for energies and data");
212 } else { 201 } else {
213 202
214 delete energies; 203 delete energies;
215 energies = dataX; 204 energies = dataX;
216 205
217 delete data; 206 delete data;
218 data = dataY; 207 data = dataY;
219 //G4cout << "Size of energies: " << ener 208 //G4cout << "Size of energies: " << energies->size() << G4endl
220 //<< "Size of data: " << data->size() << 209 //<< "Size of data: " << data->size() << G4endl;
221 } 210 }
222 } 211 }
223 } 212 }
224 213
225 //....oooOO0OOooo........oooOO0OOooo........oo <<
226 <<
227 void G4EMDataSet::SetLogEnergiesData(G4DataVec 214 void G4EMDataSet::SetLogEnergiesData(G4DataVector* dataX,
228 G4DataVec 215 G4DataVector* dataY,
229 G4DataVec 216 G4DataVector* data_logX,
230 G4DataVector* data_logY, 217 G4DataVector* data_logY,
231 G4int /* 218 G4int /* componentId */)
232 { 219 {
233 if(!dataX || !dataY || !data_logX || !data_l 220 if(!dataX || !dataY || !data_logX || !data_logY) {
234 G4Exception("G4EMDataSet::SetEnergiesData" 221 G4Exception("G4EMDataSet::SetEnergiesData",
235 "em1012",FatalException,"new interpola 222 "em1012",FatalException,"new interpolation == 0");
236 } else { 223 } else {
237 if (dataX->size() != dataY->size() || 224 if (dataX->size() != dataY->size() ||
238 dataX->size() != data_logX->size() || 225 dataX->size() != data_logX->size() ||
239 dataX->size() != data_logY->size()) { 226 dataX->size() != data_logY->size()) {
240 G4Exception("G4EMDataSet::SetEnergiesDat 227 G4Exception("G4EMDataSet::SetEnergiesData",
241 "em1012",FatalException,"different size 228 "em1012",FatalException,"different size for energies and data");
242 } else { 229 } else {
243 230
244 delete energies; 231 delete energies;
245 energies = dataX; 232 energies = dataX;
246 233
247 delete data; 234 delete data;
248 data = dataY; 235 data = dataY;
249 236
250 delete log_energies; 237 delete log_energies;
251 log_energies = data_logX; 238 log_energies = data_logX;
252 239
253 delete log_data; 240 delete log_data;
254 log_data = data_logY; 241 log_data = data_logY;
255 //G4cout << "Size of energies: " << ener 242 //G4cout << "Size of energies: " << energies->size() << G4endl
256 //<< "Size of data: " << data->size() << 243 //<< "Size of data: " << data->size() << G4endl;
257 } 244 }
258 } 245 }
259 } 246 }
260 247
261 //....oooOO0OOooo........oooOO0OOooo........oo <<
262 <<
263 G4bool G4EMDataSet::LoadData(const G4String& f 248 G4bool G4EMDataSet::LoadData(const G4String& fileName)
264 { 249 {
265 // The file is organized into four columns: 250 // The file is organized into four columns:
266 // 1st column contains the values of energy 251 // 1st column contains the values of energy
267 // 2nd column contains the corresponding data 252 // 2nd column contains the corresponding data value
268 // The file terminates with the pattern: -1 253 // The file terminates with the pattern: -1 -1
269 // -2 254 // -2 -2
270 255
271 G4String fullFileName(FullFileName(fileName) 256 G4String fullFileName(FullFileName(fileName));
272 std::ifstream in(fullFileName); 257 std::ifstream in(fullFileName);
273 258
274 if (!in.is_open()) 259 if (!in.is_open())
275 { 260 {
276 G4String message("data file \""); 261 G4String message("data file \"");
277 message += fullFileName; 262 message += fullFileName;
278 message += "\" not found"; 263 message += "\" not found";
279 G4Exception("G4EMDataSet::LoadData", 264 G4Exception("G4EMDataSet::LoadData",
280 "em1012",FatalException,message); 265 "em1012",FatalException,message);
281 return false; 266 return false;
282 } 267 }
283 268
284 delete energies; 269 delete energies;
285 delete data; 270 delete data;
286 delete log_energies; 271 delete log_energies;
287 delete log_data; 272 delete log_data;
288 energies = new G4DataVector; 273 energies = new G4DataVector;
289 data = new G4DataVector; 274 data = new G4DataVector;
290 log_energies = new G4DataVector; 275 log_energies = new G4DataVector;
291 log_data = new G4DataVector; 276 log_data = new G4DataVector;
292 277
293 G4double a, b; 278 G4double a, b;
>> 279 //G4int k = 0;
>> 280 //G4int nColumns = 2;
>> 281
294 do 282 do
295 { 283 {
296 in >> a >> b; 284 in >> a >> b;
297 285
298 if (a != -1 && a != -2) 286 if (a != -1 && a != -2)
299 { 287 {
300 if (a==0.) { a=1e-300; } 288 if (a==0.) { a=1e-300; }
301 if (b==0.) { b=1e-300; } 289 if (b==0.) { b=1e-300; }
302 a *= unitEnergies; 290 a *= unitEnergies;
303 b *= unitData; 291 b *= unitData;
304 energies->push_back(a); 292 energies->push_back(a);
305 log_energies->push_back(std::log10(a)); 293 log_energies->push_back(std::log10(a));
306 data->push_back(b); 294 data->push_back(b);
307 log_data->push_back(std::log10(b)); 295 log_data->push_back(std::log10(b));
308 } 296 }
309 } 297 }
310 while (a != -2); 298 while (a != -2);
311 299
312 if (randomSet) { BuildPdf(); } 300 if (randomSet) { BuildPdf(); }
313 301
314 return true; 302 return true;
315 } 303 }
316 304
317 //....oooOO0OOooo........oooOO0OOooo........oo <<
318 305
319 G4bool G4EMDataSet::LoadNonLogData(const G4Str 306 G4bool G4EMDataSet::LoadNonLogData(const G4String& fileName)
320 { 307 {
321 // The file is organized into four columns: 308 // The file is organized into four columns:
322 // 1st column contains the values of energy 309 // 1st column contains the values of energy
323 // 2nd column contains the corresponding data 310 // 2nd column contains the corresponding data value
324 // The file terminates with the pattern: -1 311 // The file terminates with the pattern: -1 -1
325 // -2 312 // -2 -2
326 313
327 G4String fullFileName(FullFileName(fileName) 314 G4String fullFileName(FullFileName(fileName));
328 std::ifstream in(fullFileName); 315 std::ifstream in(fullFileName);
329 if (!in.is_open()) 316 if (!in.is_open())
330 { 317 {
331 G4String message("data file \""); 318 G4String message("data file \"");
332 message += fullFileName; 319 message += fullFileName;
333 message += "\" not found"; 320 message += "\" not found";
334 G4Exception("G4EMDataSet::LoadNonLogData 321 G4Exception("G4EMDataSet::LoadNonLogData",
335 "em1012",FatalException,message); 322 "em1012",FatalException,message);
336 } 323 }
337 324
338 G4DataVector* argEnergies=new G4DataVector; 325 G4DataVector* argEnergies=new G4DataVector;
339 G4DataVector* argData=new G4DataVector; 326 G4DataVector* argData=new G4DataVector;
340 327
341 G4double a; 328 G4double a;
342 G4int k = 0; 329 G4int k = 0;
343 G4int nColumns = 2; 330 G4int nColumns = 2;
344 331
345 do 332 do
346 { 333 {
347 in >> a; 334 in >> a;
348 335
349 if (a != -1 && a != -2) 336 if (a != -1 && a != -2)
350 { 337 {
351 if (k%nColumns == 0) 338 if (k%nColumns == 0)
352 { 339 {
353 argEnergies->push_back(a*unitEnergies); 340 argEnergies->push_back(a*unitEnergies);
354 } 341 }
355 else if (k%nColumns == 1) 342 else if (k%nColumns == 1)
356 { 343 {
357 argData->push_back(a*unitData); 344 argData->push_back(a*unitData);
358 } 345 }
359 k++; 346 k++;
360 } 347 }
361 } 348 }
362 while (a != -2); 349 while (a != -2);
363 350
364 SetEnergiesData(argEnergies, argData, 0); 351 SetEnergiesData(argEnergies, argData, 0);
365 352
366 if (randomSet) BuildPdf(); 353 if (randomSet) BuildPdf();
367 354
368 return true; 355 return true;
369 } 356 }
370 357
371 //....oooOO0OOooo........oooOO0OOooo........oo << 358
372 359
373 G4bool G4EMDataSet::SaveData(const G4String& n 360 G4bool G4EMDataSet::SaveData(const G4String& name) const
374 { 361 {
375 // The file is organized into two columns: 362 // The file is organized into two columns:
376 // 1st column is the energy 363 // 1st column is the energy
377 // 2nd column is the corresponding value 364 // 2nd column is the corresponding value
378 // The file terminates with the pattern: -1 365 // The file terminates with the pattern: -1 -1
379 // -2 366 // -2 -2
380 367
381 G4String fullFileName(FullFileName(name)); 368 G4String fullFileName(FullFileName(name));
382 std::ofstream out(fullFileName); 369 std::ofstream out(fullFileName);
383 370
384 if (!out.is_open()) 371 if (!out.is_open())
385 { 372 {
386 G4String message("cannot open \""); 373 G4String message("cannot open \"");
387 message+=fullFileName; 374 message+=fullFileName;
388 message+="\""; 375 message+="\"";
389 G4Exception("G4EMDataSet::SaveData", 376 G4Exception("G4EMDataSet::SaveData",
390 "em1012",FatalException,message); 377 "em1012",FatalException,message);
391 } 378 }
392 379
393 out.precision(10); 380 out.precision(10);
394 out.width(15); 381 out.width(15);
395 out.setf(std::ofstream::left); 382 out.setf(std::ofstream::left);
396 383
397 if (energies!=0 && data!=0) 384 if (energies!=0 && data!=0)
398 { 385 {
399 G4DataVector::const_iterator i(energies- 386 G4DataVector::const_iterator i(energies->begin());
400 G4DataVector::const_iterator endI(energi 387 G4DataVector::const_iterator endI(energies->end());
401 G4DataVector::const_iterator j(data->beg 388 G4DataVector::const_iterator j(data->begin());
402 389
403 while (i!=endI) 390 while (i!=endI)
404 { 391 {
405 out.precision(10); 392 out.precision(10);
406 out.width(15); 393 out.width(15);
407 out.setf(std::ofstream::left); 394 out.setf(std::ofstream::left);
408 out << ((*i)/unitEnergies) << ' '; 395 out << ((*i)/unitEnergies) << ' ';
409 396
410 out.precision(10); 397 out.precision(10);
411 out.width(15); 398 out.width(15);
412 out.setf(std::ofstream::left); 399 out.setf(std::ofstream::left);
413 out << ((*j)/unitData) << std::endl; 400 out << ((*j)/unitData) << std::endl;
414 401
415 i++; 402 i++;
416 j++; 403 j++;
417 } 404 }
418 } 405 }
419 406
420 out.precision(10); 407 out.precision(10);
421 out.width(15); 408 out.width(15);
422 out.setf(std::ofstream::left); 409 out.setf(std::ofstream::left);
423 out << -1.f << ' '; 410 out << -1.f << ' ';
424 411
425 out.precision(10); 412 out.precision(10);
426 out.width(15); 413 out.width(15);
427 out.setf(std::ofstream::left); 414 out.setf(std::ofstream::left);
428 out << -1.f << std::endl; 415 out << -1.f << std::endl;
429 416
430 out.precision(10); 417 out.precision(10);
431 out.width(15); 418 out.width(15);
432 out.setf(std::ofstream::left); 419 out.setf(std::ofstream::left);
433 out << -2.f << ' '; 420 out << -2.f << ' ';
434 421
435 out.precision(10); 422 out.precision(10);
436 out.width(15); 423 out.width(15);
437 out.setf(std::ofstream::left); 424 out.setf(std::ofstream::left);
438 out << -2.f << std::endl; 425 out << -2.f << std::endl;
439 426
440 return true; 427 return true;
441 } 428 }
442 429
443 //....oooOO0OOooo........oooOO0OOooo........oo <<
444 430
445 std::size_t G4EMDataSet::FindLowerBound(G4doub << 431
>> 432 size_t G4EMDataSet::FindLowerBound(G4double x) const
446 { 433 {
447 std::size_t lowerBound = 0; << 434 size_t lowerBound = 0;
448 std::size_t upperBound(energies->size() - 1) << 435 size_t upperBound(energies->size() - 1);
449 436
450 while (lowerBound <= upperBound) 437 while (lowerBound <= upperBound)
451 { 438 {
452 std::size_t midBin((lowerBound + upperBo << 439 size_t midBin((lowerBound + upperBound) / 2);
453 440
454 if (x < (*energies)[midBin]) upperBound 441 if (x < (*energies)[midBin]) upperBound = midBin - 1;
455 else lowerBound = midBin + 1; 442 else lowerBound = midBin + 1;
456 } 443 }
457 444
458 return upperBound; 445 return upperBound;
459 } 446 }
460 447
461 //....oooOO0OOooo........oooOO0OOooo........oo <<
462 448
463 std::size_t G4EMDataSet::FindLowerBound(G4doub << 449 size_t G4EMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
464 { 450 {
465 std::size_t lowerBound = 0;; << 451 size_t lowerBound = 0;;
466 std::size_t upperBound(values->size() - 1); << 452 size_t upperBound(values->size() - 1);
467 453
468 while (lowerBound <= upperBound) 454 while (lowerBound <= upperBound)
469 { 455 {
470 std::size_t midBin((lowerBound + upperBo << 456 size_t midBin((lowerBound + upperBound) / 2);
471 457
472 if (x < (*values)[midBin]) upperBound = 458 if (x < (*values)[midBin]) upperBound = midBin - 1;
473 else lowerBound = midBin + 1; 459 else lowerBound = midBin + 1;
474 } 460 }
475 461
476 return upperBound; 462 return upperBound;
477 } 463 }
478 464
479 //....oooOO0OOooo........oooOO0OOooo........oo <<
480 465
481 G4String G4EMDataSet::FullFileName(const G4Str 466 G4String G4EMDataSet::FullFileName(const G4String& name) const
482 { 467 {
483 const char* path = G4FindDataDir("G4LEDATA") << 468 char* path = getenv("G4LEDATA");
484 if (!path) { 469 if (!path) {
485 G4Exception("G4EMDataSet::FullFileName", 470 G4Exception("G4EMDataSet::FullFileName",
486 "em0006",FatalException,"G4LEDATA envi 471 "em0006",FatalException,"G4LEDATA environment variable not set");
487 return ""; 472 return "";
488 } 473 }
489 std::ostringstream fullFileName; 474 std::ostringstream fullFileName;
490 fullFileName << path << '/' << name << z << 475 fullFileName << path << '/' << name << z << ".dat";
491 476
492 return G4String(fullFileName.str().c_str()); 477 return G4String(fullFileName.str().c_str());
493 } 478 }
494 479
495 //....oooOO0OOooo........oooOO0OOooo........oo << 480
496 481
497 void G4EMDataSet::BuildPdf() 482 void G4EMDataSet::BuildPdf()
498 { 483 {
499 pdf = new G4DataVector; 484 pdf = new G4DataVector;
500 G4Integrator <G4EMDataSet, G4double(G4EMData 485 G4Integrator <G4EMDataSet, G4double(G4EMDataSet::*)(G4double)> integrator;
501 486
502 std::size_t nData = data->size(); << 487 G4int nData = data->size();
503 pdf->push_back(0.); 488 pdf->push_back(0.);
504 489
505 // Integrate the data distribution 490 // Integrate the data distribution
506 std::size_t i; << 491 G4int i;
507 G4double totalSum = 0.; 492 G4double totalSum = 0.;
508 for (i=1; i<nData; ++i) << 493 for (i=1; i<nData; i++)
509 { 494 {
510 G4double xLow = (*energies)[i-1]; 495 G4double xLow = (*energies)[i-1];
511 G4double xHigh = (*energies)[i]; 496 G4double xHigh = (*energies)[i];
512 G4double sum = integrator.Legendre96(thi 497 G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh);
513 totalSum = totalSum + sum; 498 totalSum = totalSum + sum;
514 pdf->push_back(totalSum); 499 pdf->push_back(totalSum);
515 } 500 }
516 501
517 // Normalize to the last bin 502 // Normalize to the last bin
518 G4double tot = 0.; 503 G4double tot = 0.;
519 if (totalSum > 0.) tot = 1. / totalSum; 504 if (totalSum > 0.) tot = 1. / totalSum;
520 for (i=1; i<nData; ++i) << 505 for (i=1; i<nData; i++)
521 { 506 {
522 (*pdf)[i] = (*pdf)[i] * tot; 507 (*pdf)[i] = (*pdf)[i] * tot;
523 } 508 }
524 } 509 }
525 510
526 //....oooOO0OOooo........oooOO0OOooo........oo <<
527 <<
528 G4double G4EMDataSet::RandomSelect(G4int /* co 511 G4double G4EMDataSet::RandomSelect(G4int /* componentId */) const
529 { 512 {
530 G4double value = 0.; 513 G4double value = 0.;
531 // Random select a X value according to the 514 // Random select a X value according to the cumulative probability distribution
532 // derived from the data 515 // derived from the data
533 516
534 if (!pdf) { 517 if (!pdf) {
535 G4Exception("G4EMDataSet::RandomSelect", 518 G4Exception("G4EMDataSet::RandomSelect",
536 "em1012",FatalException,"PDF has not b 519 "em1012",FatalException,"PDF has not been created for this data set");
537 return value; 520 return value;
538 } 521 }
539 522
540 G4double x = G4UniformRand(); 523 G4double x = G4UniformRand();
541 524
542 // Locate the random value in the X vector b 525 // Locate the random value in the X vector based on the PDF
543 G4int bin = (G4int)FindLowerBound(x,pdf); << 526 size_t bin = FindLowerBound(x,pdf);
544 527
545 // Interpolate the PDF to calculate the X va 528 // Interpolate the PDF to calculate the X value:
546 // linear interpolation in the first bin (to 529 // linear interpolation in the first bin (to avoid problem with 0),
547 // interpolation with associated data set al 530 // interpolation with associated data set algorithm in other bins
>> 531
548 G4LinInterpolation linearAlgo; 532 G4LinInterpolation linearAlgo;
549 if (bin == 0) value = linearAlgo.Calculate(x 533 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
550 else value = algorithm->Calculate(x, bin, *p 534 else value = algorithm->Calculate(x, bin, *pdf, *energies);
551 535
>> 536 // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
552 return value; 537 return value;
553 } 538 }
554 539
555 //....oooOO0OOooo........oooOO0OOooo........oo <<
556 <<
557 G4double G4EMDataSet::IntegrationFunction(G4do 540 G4double G4EMDataSet::IntegrationFunction(G4double x)
558 { 541 {
559 // This function is needed by G4Integrator t 542 // This function is needed by G4Integrator to calculate the integral of the data distribution
>> 543
560 G4double y = 0; 544 G4double y = 0;
561 545
562 // Locate the random value in the X vector b << 546 // Locate the random value in the X vector based on the PDF
563 G4int bin = (G4int)FindLowerBound(x); << 547 size_t bin = FindLowerBound(x);
564 548
565 // Interpolate to calculate the X value: 549 // Interpolate to calculate the X value:
566 // linear interpolation in the first bin (to 550 // linear interpolation in the first bin (to avoid problem with 0),
567 // interpolation with associated algorithm i 551 // interpolation with associated algorithm in other bins
568 552
569 G4LinInterpolation linearAlgo; 553 G4LinInterpolation linearAlgo;
570 554
571 if (bin == 0) y = linearAlgo.Calculate(x, bi 555 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
572 else y = algorithm->Calculate(x, bin, *energ 556 else y = algorithm->Calculate(x, bin, *energies, *data);
573 557
574 return y; 558 return y;
575 } 559 }
576 560
577 561