| Geant4 Cross Reference |
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
45 G4RDEMDataSet::G4RDEMDataSet(G4int Z,
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
65 G4RDEMDataSet::G4RDEMDataSet(G4int argZ,
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
98 G4RDEMDataSet::~G4RDEMDataSet()
99 {
100 delete algorithm;
101 if (energies) delete energies;
102 if (data) delete data;
103 if (pdf) delete pdf;
104 }
105
106 G4double G4RDEMDataSet::FindValue(G4double energy, G4int /* componentId */) const
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
141 void G4RDEMDataSet::SetEnergiesData(G4DataVector* dataX,
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
162 G4bool G4RDEMDataSet::LoadData(const G4String& fileName)
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
209 G4bool G4RDEMDataSet::SaveData(const G4String& name) const
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
279 size_t G4RDEMDataSet::FindLowerBound(G4double x) const
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
296 size_t G4RDEMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
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
313 G4String G4RDEMDataSet::FullFileName(const G4String& name) const
314 {
315 const char* path = G4FindDataDir("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
327 void G4RDEMDataSet::BuildPdf()
328 {
329 pdf = new G4DataVector;
330 G4Integrator <G4RDEMDataSet, G4double(G4RDEMDataSet::*)(G4double)> integrator;
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
357 G4double G4RDEMDataSet::RandomSelect(G4int /* componentId */) const
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.;
366 G4double x = G4UniformRand();
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
383 G4double G4RDEMDataSet::IntegrationFunction(G4double x)
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
404