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