ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4EMDataSet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4EMDataSet.cc
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;
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 
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 
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 
467 {
468  char* path = std::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 
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 
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 
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