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