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