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