source: trunk/source/processes/electromagnetic/lowenergy/src/G4EMDataSet.cc@ 1229

Last change on this file since 1229 was 1228, checked in by garnier, 16 years ago

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