ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4DataSet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4DataSet.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 //
27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 //
30 // History:
31 // -----------
32 // 31 Jul 2001 MGP Created
33 // 31 Jul 2008 MGP Revised and renamed to G4DataSet
34 //
35 // -------------------------------------------------------------------
36 
37 #include "G4DataSet.hh"
38 #include "G4IInterpolator.hh"
39 #include <fstream>
40 #include <sstream>
41 #include "G4Integrator.hh"
42 #include "Randomize.hh"
43 #include "G4LinInterpolation.hh"
44 
45 
47  G4IInterpolator* algo,
48  G4double xUnit,
49  G4double yUnit,
50  G4bool random):
51  z(Z),
52  energies(0),
53  data(0),
54  algorithm(algo),
55  unitEnergies(xUnit),
56  unitData(yUnit),
57  pdf(0),
58  randomSet(random)
59 {
60  if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
61  "pii00000101",
63  "Interpolation == 0");
64  if (randomSet) BuildPdf();
65 }
66 
68  G4DataVector* dataX,
69  G4DataVector* dataY,
70  G4IInterpolator* algo,
71  G4double xUnit,
72  G4double yUnit,
73  G4bool random):
74  z(argZ),
75  energies(dataX),
76  data(dataY),
77  algorithm(algo),
78  unitEnergies(xUnit),
79  unitData(yUnit),
80  pdf(0),
81  randomSet(random)
82 {
83  if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
84  "pii00000110",
86  "Interpolation == 0");
87 
88  if ((energies == 0) ^ (data == 0))
89  G4Exception("G4DataSet::G4DataSet",
90  "pii00000111-",
92  "different size for energies and data (zero case)");
93 
94  if (energies == 0) return;
95 
96  if (energies->size() != data->size())
97  G4Exception("G4DataSet::G4DataSet",
98  "pii00000112",
100  "different size for energies and data");
101 
102  if (randomSet) BuildPdf();
103 }
104 
106 {
107  delete algorithm;
108  if (energies) delete energies;
109  if (data) delete data;
110  if (pdf) delete pdf;
111 }
112 
114 {
115  if (!energies) G4Exception("G4DataSet::FindValue",
116  "pii00000120",
118  "energies == 0");
119  if (energies->empty()) return 0;
120  if (energy <= (*energies)[0]) return (*data)[0];
121 
122  size_t i = energies->size()-1;
123  if (energy >= (*energies)[i]) return (*data)[i];
124 
125  G4double interpolated = algorithm->Calculate(energy,FindLowerBound(energy),*energies,*data);
126  return interpolated;
127 }
128 
129 
130 void G4DataSet::PrintData(void) const
131 {
132  if (!energies)
133  {
134  G4cout << "Data not available." << G4endl;
135  }
136  else
137  {
138  size_t size = energies->size();
139  for (size_t i(0); i<size; i++)
140  {
141  G4cout << "Point: " << ((*energies)[i]/unitEnergies)
142  << " - Data value: " << ((*data)[i]/unitData);
143  if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
144  G4cout << G4endl;
145  }
146  }
147 }
148 
149 
151  G4DataVector* dataY,
152  G4int /* componentId */)
153 {
154  if (energies) delete energies;
155  energies = dataX;
156 
157  if (data) delete data;
158  data = dataY;
159 
160  if ((energies == 0) ^ (data==0))
161  G4Exception("G4DataSet::SetEnergiesData",
162  "pii00000130",
164  "different size for energies and data (zero case)");
165 
166  if (energies == 0) return;
167 
168  if (energies->size() != data->size())
169  G4Exception("G4DataSet::SetEnergiesData",
170  "pii00000131",
172  "different size for energies and data");
173 }
174 
176 {
177  // The file is organized into two columns:
178  // 1st column is the energy
179  // 2nd column is the corresponding value
180  // The file terminates with the pattern: -1 -1
181  // -2 -2
182 
183  G4String fullFileName(FullFileName(fileName));
184  std::ifstream in(fullFileName);
185 
186  if (!in.is_open())
187  {
188 
189  std::ostringstream message;
190  message << "G4DataSet::LoadData - data file " << fullFileName << " not found";
191 
192  G4Exception("G4CompositeDataSet::LoadData",
193  "pii00000140",
195  message.str().c_str());
196  }
197 
198  G4DataVector* argEnergies=new G4DataVector;
199  G4DataVector* argData=new G4DataVector;
200 
201  G4double a;
202  bool energyColumn(true);
203 
204  do
205  {
206  in >> a;
207 
208  if (a!=-1 && a!=-2)
209  {
210  if (energyColumn)
211  {
212  // std::cout << fullFileName << ", a = " << a <<std::endl;
213  argEnergies->push_back(a*unitEnergies);
214  }
215  else
216  argData->push_back(a*unitData);
217  energyColumn=(!energyColumn);
218  }
219  }
220  while (a != -2);
221 
222  SetEnergiesData(argEnergies, argData, 0);
223  if (randomSet) BuildPdf();
224 
225  return true;
226 }
227 
229 {
230  // The file is organized into two columns:
231  // 1st column is the energy
232  // 2nd column is the corresponding value
233  // The file terminates with the pattern: -1 -1
234  // -2 -2
235 
236  G4String fullFileName(FullFileName(name));
237  std::ofstream out(fullFileName);
238 
239  if (!out.is_open())
240  {
241 
242  std::ostringstream message;
243  message << "G4DataSet:: SaveData - cannot open " << fullFileName;
244 
245  G4Exception("G4CompositeDataSet::SaveData",
246  "pii00000150",
248  message.str().c_str());
249 
250  }
251 
252  out.precision(10);
253  out.width(15);
254  out.setf(std::ofstream::left);
255 
256  if (energies!=0 && data!=0)
257  {
258  G4DataVector::const_iterator i(energies->begin());
259  G4DataVector::const_iterator endI(energies->end());
260  G4DataVector::const_iterator j(data->begin());
261 
262  while (i!=endI)
263  {
264  out.precision(10);
265  out.width(15);
266  out.setf(std::ofstream::left);
267  out << ((*i)/unitEnergies) << ' ';
268 
269  out.precision(10);
270  out.width(15);
271  out.setf(std::ofstream::left);
272  out << ((*j)/unitData) << std::endl;
273 
274  i++;
275  j++;
276  }
277  }
278 
279  out.precision(10);
280  out.width(15);
281  out.setf(std::ofstream::left);
282  out << -1.f << ' ';
283 
284  out.precision(10);
285  out.width(15);
286  out.setf(std::ofstream::left);
287  out << -1.f << std::endl;
288 
289  out.precision(10);
290  out.width(15);
291  out.setf(std::ofstream::left);
292  out << -2.f << ' ';
293 
294  out.precision(10);
295  out.width(15);
296  out.setf(std::ofstream::left);
297  out << -2.f << std::endl;
298 
299  return true;
300 }
301 
303 {
304  size_t lowerBound = 0;
305  size_t upperBound(energies->size() - 1);
306 
307  while (lowerBound <= upperBound)
308  {
309  size_t midBin((lowerBound + upperBound) / 2);
310 
311  if (x < (*energies)[midBin]) upperBound = midBin - 1;
312  else lowerBound = midBin + 1;
313  }
314 
315  return upperBound;
316 }
317 
318 
320 {
321  size_t lowerBound = 0;;
322  size_t upperBound(values->size() - 1);
323 
324  while (lowerBound <= upperBound)
325  {
326  size_t midBin((lowerBound + upperBound) / 2);
327 
328  if (x < (*values)[midBin]) upperBound = midBin - 1;
329  else lowerBound = midBin + 1;
330  }
331 
332  return upperBound;
333 }
334 
335 
337 {
338  char* path = std::getenv("G4PIIDATA");
339  if (!path)
340  G4Exception("G4DataSet::FullFileName",
341  "pii00000160",
343  "G4PIIDATA environment variable not set");
344 
345  std::ostringstream fullFileName;
346  fullFileName << path << '/' << name << z << ".dat";
347 
348  return G4String(fullFileName.str().c_str());
349 }
350 
351 
353 {
354  pdf = new G4DataVector;
356 
357  G4int nData = data->size();
358  pdf->push_back(0.);
359 
360  // Integrate the data distribution
361  G4int i;
362  G4double totalSum = 0.;
363  for (i=1; i<nData; i++)
364  {
365  G4double xLow = (*energies)[i-1];
366  G4double xHigh = (*energies)[i];
367  G4double sum = integrator.Legendre96(this, &G4DataSet::IntegrationFunction, xLow, xHigh);
368  totalSum = totalSum + sum;
369  pdf->push_back(totalSum);
370  }
371 
372  // Normalize to the last bin
373  G4double tot = 0.;
374  if (totalSum > 0.) tot = 1. / totalSum;
375  for (i=1; i<nData; i++)
376  {
377  (*pdf)[i] = (*pdf)[i] * tot;
378  }
379 }
380 
381 
382 G4double G4DataSet::RandomSelect(G4int /* componentId */) const
383 {
384  // Random select a X value according to the cumulative probability distribution
385  // derived from the data
386 
387  if (!pdf) G4Exception("G4DataSet::RandomSelect",
388  "pii00000170",
390  "PDF has not been created for this data set");
391 
392  G4double value = 0.;
394 
395  // Locate the random value in the X vector based on the PDF
396  size_t bin = FindLowerBound(x,pdf);
397 
398  // Interpolate the PDF to calculate the X value:
399  // linear interpolation in the first bin (to avoid problem with 0),
400  // interpolation with associated data set algorithm in other bins
401 
402  G4LinInterpolation linearAlgo;
403  if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
404  else value = algorithm->Calculate(x, bin, *pdf, *energies);
405 
406  // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
407  return value;
408 }
409 
411 {
412  // This function is needed by G4Integrator to calculate the integral of the data distribution
413 
414  G4double y = 0;
415 
416  // Locate the random value in the X vector based on the PDF
417  size_t bin = FindLowerBound(x);
418 
419  // Interpolate to calculate the X value:
420  // linear interpolation in the first bin (to avoid problem with 0),
421  // interpolation with associated algorithm in other bins
422 
423  G4LinInterpolation linearAlgo;
424 
425  if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
426  else y = algorithm->Calculate(x, bin, *energies, *data);
427 
428  return y;
429 }
430