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