Geant4  10.01.p01
G4EMDataSet.cc
Go to the documentation of this file.
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 
63  G4VDataSetAlgorithm* algo,
64  G4double xUnit,
65  G4double yUnit,
66  G4bool random):
67  z(Z),
68  energies(0),
69  data(0),
70  log_energies(0),
71  log_data(0),
72  algorithm(algo),
73  unitEnergies(xUnit),
74  unitData(yUnit),
75  pdf(0),
76  randomSet(random)
77 {
78  if (algorithm == 0) {
79  G4Exception("G4EMDataSet::G4EMDataSet",
80  "em1012",FatalException,"interpolation == 0");
81  } else if (randomSet) { BuildPdf(); }
82 }
83 
85  G4DataVector* dataX,
86  G4DataVector* dataY,
87  G4VDataSetAlgorithm* algo,
88  G4double xUnit,
89  G4double yUnit,
90  G4bool random):
91  z(argZ),
92  energies(dataX),
93  data(dataY),
94  log_energies(0),
95  log_data(0),
96  algorithm(algo),
97  unitEnergies(xUnit),
98  unitData(yUnit),
99  pdf(0),
100  randomSet(random)
101 {
102  if (algorithm == 0) {
103  G4Exception("G4EMDataSet::G4EMDataSet",
104  "em1012",FatalException,"interpolation == 0");
105  } else if ((energies == 0) ^ (data == 0)) {
106  G4Exception("G4EMDataSet::G4EMDataSet",
107  "em1012",FatalException,"different size for energies and data (zero case)");
108  } else if (energies->size() != data->size()) {
109  G4Exception("G4EMDataSet::G4EMDataSet",
110  "em1012",FatalException,"different size for energies and data");
111  } else if (randomSet) {
112  BuildPdf();
113  }
114 }
115 
117  G4DataVector* dataX,
118  G4DataVector* dataY,
119  G4DataVector* dataLogX,
120  G4DataVector* dataLogY,
121  G4VDataSetAlgorithm* algo,
122  G4double xUnit,
123  G4double yUnit,
124  G4bool random):
125  z(argZ),
126  energies(dataX),
127  data(dataY),
128  log_energies(dataLogX),
129  log_data(dataLogY),
130  algorithm(algo),
131  unitEnergies(xUnit),
132  unitData(yUnit),
133  pdf(0),
134  randomSet(random)
135 {
136  if (algorithm == 0) {
137  G4Exception("G4EMDataSet::G4EMDataSet",
138  "em1012",FatalException,"interpolation == 0");
139  } else if ((energies == 0) ^ (data == 0)) {
140  G4Exception("G4EMDataSet::G4EMDataSet",
141  "em1012",FatalException,"different size for energies and data (zero case)");
142  } else if (energies->size() != data->size()) {
143  G4Exception("G4EMDataSet::G4EMDataSet",
144  "em1012",FatalException,"different size for energies and data");
145  } else if ((log_energies == 0) ^ (log_data == 0)) {
146  G4Exception("G4EMDataSet::G4EMDataSet",
147  "em1012",FatalException,"different size for log energies and log data (zero case)");
148  } else if (log_energies->size() != log_data->size()) {
149  G4Exception("G4EMDataSet::G4EMDataSet",
150  "em1012",FatalException,"different size for log energies and log data");
151  } else if (randomSet) {
152  BuildPdf();
153  }
154 }
155 
156 
158 {
159  delete algorithm;
160  delete energies;
161  delete data;
162  delete pdf;
163  delete log_energies;
164  delete log_data;
165 }
166 
167 
169 {
170  if (!energies) {
171  G4Exception("G4EMDataSet::FindValue",
172  "em1012",FatalException,"energies == 0");
173  return 0.0;
174  } else if (energies->empty()) {
175  return 0.0;
176  } else if (energy <= (*energies)[0]) {
177  return (*data)[0];
178  }
179  size_t i = energies->size()-1;
180  if (energy >= (*energies)[i]) { return (*data)[i]; }
181 
182  //Nicolas A. Karakatsanis: Check if the logarithmic data have been loaded to decide
183  // which Interpolation-Calculation method will be applied
184  if (log_energies != 0)
185  {
186  return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data, *log_energies, *log_data);
187  }
188 
189  return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
190 }
191 
192 
193 void G4EMDataSet::PrintData(void) const
194 {
195  if (!energies)
196  {
197  G4cout << "Data not available." << G4endl;
198  }
199  else
200  {
201  size_t size = energies->size();
202  for (size_t i(0); i<size; i++)
203  {
204  G4cout << "Point: " << ((*energies)[i]/unitEnergies)
205  << " - Data value: " << ((*data)[i]/unitData);
206  if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
207  G4cout << G4endl;
208  }
209  }
210 }
211 
212 
214  G4DataVector* dataY,
215  G4int /* componentId */)
216 {
217  if (energies) { delete energies; }
218  energies = dataX;
219 
220  if (data) { delete data; }
221  data = dataY;
222 
223  if ((energies == 0) ^ (data==0)) {
224  G4Exception("G4EMDataSet::SetEnergiesData",
225  "em1012",FatalException,"different size for energies and data (zero case)");
226  return;
227  } else if (energies == 0) { return; }
228 
229  //G4cout << "Size of energies: " << energies->size() << G4endl << "Size of data: " << data->size() << G4endl;
230  if (energies->size() != data->size()) {
231  G4Exception("G4EMDataSet::SetEnergiesData",
232  "em1012",FatalException,"different size for energies and data");
233  }
234 }
235 
237  G4DataVector* dataY,
238  G4DataVector* data_logX,
239  G4DataVector* data_logY,
240  G4int /* componentId */)
241 {
242  //Load of the actual energy and data values
243  if (energies) { delete energies; }
244  energies = dataX;
245  if (data) { delete data; }
246  data = dataY;
247  //Load of the logarithmic energy and data values
248  if (log_energies) { delete log_energies; }
249  log_energies = data_logX;
250  if (log_data) { delete log_data; }
251  log_data = data_logY;
252 
253  //Check if data loaded properly from data files
254  if ( !energies ) {
255  if(data || log_energies || log_data ) {
256  G4Exception("G4EMDataSet::SetLogEnergiesData",
257  "em1012",FatalException,"inconsistent data");
258  }
259  return;
260  } else {
261  if ( !data ) {
262  G4Exception("G4EMDataSet::SetLogEnergiesData",
263  "em1012",FatalException,"only energy, no data");
264  return;
265  } else if (energies->size() != data->size()) {
266  G4Exception("G4EMDataSet::SetLogEnergiesData",
267  "em1012",FatalException,"different size for energies and data");
268  return;
269  }
270  //G4cout << "Size of energies: " << energies->size() << G4endl << "Size of data: " << data->size() << G4endl << G4endl;
271 
272  //Check if logarithmic data loaded properly from data files
273  if ( !log_energies ) {
274  if(log_data) {
275  G4Exception("G4EMDataSet::SetLogEnergiesData",
276  "em1012",FatalException,"inconsistence of log_data");
277  }
278  return;
279  } else {
280  if ( !log_data ) {
281  G4Exception("G4EMDataSet::SetLogEnergiesData",
282  "em1012",FatalException,"only log_energies, no data");
283  } else if ((log_energies->size() != log_data->size()) || (log_energies->size() != data->size())) {
284  G4Exception("G4EMDataSet::SetLogEnergiesData",
285  "em1012",FatalException,"different size for log energies and data");
286  }
287  }
288  }
289  //G4cout << "Size of log energies: " << log_energies->size() << G4endl << "Size of log data: " << log_data->size() << G4endl;
290 }
291 
293 {
294  // The file is organized into four columns:
295  // 1st column contains the values of energy
296  // 2nd column contains the corresponding data value
297  // The file terminates with the pattern: -1 -1
298  // -2 -2
299 
300  G4String fullFileName(FullFileName(fileName));
301  std::ifstream in(fullFileName);
302 
303  if (!in.is_open())
304  {
305  G4String message("data file \"");
306  message += fullFileName;
307  message += "\" not found";
308  G4Exception("G4EMDataSet::LoadData",
309  "em1012",FatalException,message);
310  return false;
311  }
312 
313  delete energies;
314  delete data;
315  delete log_energies;
316  delete log_data;
317  energies = new G4DataVector;
318  data = new G4DataVector;
320  log_data = new G4DataVector;
321 
322  G4double a, b;
323  //G4int k = 0;
324  //G4int nColumns = 2;
325 
326  do
327  {
328  in >> a >> b;
329 
330  if (a != -1 && a != -2)
331  {
332  if (a==0.) { a=1e-300; }
333  if (b==0.) { b=1e-300; }
334  a *= unitEnergies;
335  b *= unitData;
336  energies->push_back(a);
337  log_energies->push_back(std::log10(a));
338  data->push_back(b);
339  log_data->push_back(std::log10(b));
340  }
341  }
342  while (a != -2);
343 
344  if (randomSet) { BuildPdf(); }
345 
346  return true;
347 }
348 
349 
351 {
352  // The file is organized into four columns:
353  // 1st column contains the values of energy
354  // 2nd column contains the corresponding data value
355  // The file terminates with the pattern: -1 -1
356  // -2 -2
357 
358  G4String fullFileName(FullFileName(fileName));
359  std::ifstream in(fullFileName);
360  if (!in.is_open())
361  {
362  G4String message("data file \"");
363  message += fullFileName;
364  message += "\" not found";
365  G4Exception("G4EMDataSet::LoadNonLogData",
366  "em1012",FatalException,message);
367  }
368 
369  G4DataVector* argEnergies=new G4DataVector;
370  G4DataVector* argData=new G4DataVector;
371 
372  G4double a;
373  G4int k = 0;
374  G4int nColumns = 2;
375 
376  do
377  {
378  in >> a;
379 
380  if (a != -1 && a != -2)
381  {
382  if (k%nColumns == 0)
383  {
384  argEnergies->push_back(a*unitEnergies);
385  }
386  else if (k%nColumns == 1)
387  {
388  argData->push_back(a*unitData);
389  }
390  k++;
391  }
392  }
393  while (a != -2);
394 
395  SetEnergiesData(argEnergies, argData, 0);
396 
397  if (randomSet) BuildPdf();
398 
399  return true;
400 }
401 
402 
403 
405 {
406  // The file is organized into two columns:
407  // 1st column is the energy
408  // 2nd column is the corresponding value
409  // The file terminates with the pattern: -1 -1
410  // -2 -2
411 
412  G4String fullFileName(FullFileName(name));
413  std::ofstream out(fullFileName);
414 
415  if (!out.is_open())
416  {
417  G4String message("cannot open \"");
418  message+=fullFileName;
419  message+="\"";
420  G4Exception("G4EMDataSet::SaveData",
421  "em1012",FatalException,message);
422  }
423 
424  out.precision(10);
425  out.width(15);
426  out.setf(std::ofstream::left);
427 
428  if (energies!=0 && data!=0)
429  {
430  G4DataVector::const_iterator i(energies->begin());
431  G4DataVector::const_iterator endI(energies->end());
432  G4DataVector::const_iterator j(data->begin());
433 
434  while (i!=endI)
435  {
436  out.precision(10);
437  out.width(15);
438  out.setf(std::ofstream::left);
439  out << ((*i)/unitEnergies) << ' ';
440 
441  out.precision(10);
442  out.width(15);
443  out.setf(std::ofstream::left);
444  out << ((*j)/unitData) << std::endl;
445 
446  i++;
447  j++;
448  }
449  }
450 
451  out.precision(10);
452  out.width(15);
453  out.setf(std::ofstream::left);
454  out << -1.f << ' ';
455 
456  out.precision(10);
457  out.width(15);
458  out.setf(std::ofstream::left);
459  out << -1.f << std::endl;
460 
461  out.precision(10);
462  out.width(15);
463  out.setf(std::ofstream::left);
464  out << -2.f << ' ';
465 
466  out.precision(10);
467  out.width(15);
468  out.setf(std::ofstream::left);
469  out << -2.f << std::endl;
470 
471  return true;
472 }
473 
474 
475 
477 {
478  size_t lowerBound = 0;
479  size_t upperBound(energies->size() - 1);
480 
481  while (lowerBound <= upperBound)
482  {
483  size_t midBin((lowerBound + upperBound) / 2);
484 
485  if (x < (*energies)[midBin]) upperBound = midBin - 1;
486  else lowerBound = midBin + 1;
487  }
488 
489  return upperBound;
490 }
491 
492 
494 {
495  size_t lowerBound = 0;;
496  size_t upperBound(values->size() - 1);
497 
498  while (lowerBound <= upperBound)
499  {
500  size_t midBin((lowerBound + upperBound) / 2);
501 
502  if (x < (*values)[midBin]) upperBound = midBin - 1;
503  else lowerBound = midBin + 1;
504  }
505 
506  return upperBound;
507 }
508 
509 
511 {
512  char* path = getenv("G4LEDATA");
513  if (!path) {
514  G4Exception("G4EMDataSet::FullFileName",
515  "em0006",FatalException,"G4LEDATA environment variable not set");
516  return "";
517  }
518  std::ostringstream fullFileName;
519  fullFileName << path << '/' << name << z << ".dat";
520 
521  return G4String(fullFileName.str().c_str());
522 }
523 
524 
525 
527 {
528  pdf = new G4DataVector;
530 
531  G4int nData = data->size();
532  pdf->push_back(0.);
533 
534  // Integrate the data distribution
535  G4int i;
536  G4double totalSum = 0.;
537  for (i=1; i<nData; i++)
538  {
539  G4double xLow = (*energies)[i-1];
540  G4double xHigh = (*energies)[i];
541  G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh);
542  totalSum = totalSum + sum;
543  pdf->push_back(totalSum);
544  }
545 
546  // Normalize to the last bin
547  G4double tot = 0.;
548  if (totalSum > 0.) tot = 1. / totalSum;
549  for (i=1; i<nData; i++)
550  {
551  (*pdf)[i] = (*pdf)[i] * tot;
552  }
553 }
554 
555 G4double G4EMDataSet::RandomSelect(G4int /* componentId */) const
556 {
557  G4double value = 0.;
558  // Random select a X value according to the cumulative probability distribution
559  // derived from the data
560 
561  if (!pdf) {
562  G4Exception("G4EMDataSet::RandomSelect",
563  "em1012",FatalException,"PDF has not been created for this data set");
564  return value;
565  }
566 
567  G4double x = G4UniformRand();
568 
569  // Locate the random value in the X vector based on the PDF
570  size_t bin = FindLowerBound(x,pdf);
571 
572  // Interpolate the PDF to calculate the X value:
573  // linear interpolation in the first bin (to avoid problem with 0),
574  // interpolation with associated data set algorithm in other bins
575 
576  G4LinInterpolation linearAlgo;
577  if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
578  else value = algorithm->Calculate(x, bin, *pdf, *energies);
579 
580  // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
581  return value;
582 }
583 
585 {
586  // This function is needed by G4Integrator to calculate the integral of the data distribution
587 
588  G4double y = 0;
589 
590  // Locate the random value in the X vector based on the PDF
591  size_t bin = FindLowerBound(x);
592 
593  // Interpolate to calculate the X value:
594  // linear interpolation in the first bin (to avoid problem with 0),
595  // interpolation with associated algorithm in other bins
596 
597  G4LinInterpolation linearAlgo;
598 
599  if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
600  else y = algorithm->Calculate(x, bin, *energies, *data);
601 
602  return y;
603 }
604 
G4DataVector * data
Definition: G4EMDataSet.hh:137
G4double Legendre96(T &typeT, F f, G4double a, G4double b)
G4double Calculate(G4double point, G4int bin, const G4DataVector &energies, const G4DataVector &data) const
G4double z
Definition: TRTMaterials.hh:39
G4String name
Definition: TRTMaterials.hh:40
virtual void SetEnergiesData(G4DataVector *xData, G4DataVector *data, G4int componentId)
Definition: G4EMDataSet.cc:213
G4DataVector * pdf
Definition: G4EMDataSet.hh:146
G4double a
Definition: TRTMaterials.hh:39
G4DataVector * log_data
Definition: G4EMDataSet.hh:139
int G4int
Definition: G4Types.hh:78
virtual G4double Calculate(G4double point, G4int bin, const G4DataVector &energies, const G4DataVector &data) const =0
virtual ~G4EMDataSet()
Definition: G4EMDataSet.cc:157
G4double unitEnergies
Definition: G4EMDataSet.hh:143
G4DataVector * energies
Definition: G4EMDataSet.hh:136
size_t FindLowerBound(G4double energy) const
Definition: G4EMDataSet.cc:476
G4double IntegrationFunction(G4double x)
Definition: G4EMDataSet.cc:584
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
G4DataVector * log_energies
Definition: G4EMDataSet.hh:138
virtual G4double FindValue(G4double x, G4int componentId=0) const
Definition: G4EMDataSet.cc:168
virtual G4bool SaveData(const G4String &fileName) const
Definition: G4EMDataSet.cc:404
G4bool randomSet
Definition: G4EMDataSet.hh:147
G4String FullFileName(const G4String &fileName) const
Definition: G4EMDataSet.cc:510
virtual void PrintData(void) const
Definition: G4EMDataSet.cc:193
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4VDataSetAlgorithm * algorithm
Definition: G4EMDataSet.hh:141
virtual G4bool LoadData(const G4String &fileName)
Definition: G4EMDataSet.cc:292
G4double energy(const ThreeVector &p, const G4double m)
virtual G4double RandomSelect(G4int componentId=0) const
Definition: G4EMDataSet.cc:555
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double unitData
Definition: G4EMDataSet.hh:144
virtual void SetLogEnergiesData(G4DataVector *xData, G4DataVector *data, G4DataVector *xLogData, G4DataVector *Logdata, G4int componentId)
Definition: G4EMDataSet.cc:236
virtual void BuildPdf()
Definition: G4EMDataSet.cc:526
virtual G4bool LoadNonLogData(const G4String &fileName)
Definition: G4EMDataSet.cc:350