| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@
29 //
30 // History:
31 // -----------
32 // 31 Jul 2001 MGP Created
33 // 31 Jul 2008 MGP Revised and rename
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
46 G4DataSet::G4DataSet(G4int Z,
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::
61 "pii00000101",
62 FatalException,
63 "Interpolation == 0");
64 if (randomSet) BuildPdf();
65 }
66
67 G4DataSet::G4DataSet(G4int argZ,
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::
84 "pii00000110",
85 FatalException,
86 "Interpolation == 0");
87
88 if ((energies == 0) ^ (data == 0))
89 G4Exception("G4DataSet::G4DataSet",
90 "pii00000111-",
91 FatalException,
92 "different size for energies and data (zer
93
94 if (energies == 0) return;
95
96 if (energies->size() != data->size())
97 G4Exception("G4DataSet::G4DataSet",
98 "pii00000112",
99 FatalException,
100 "different size for energies and data");
101
102 if (randomSet) BuildPdf();
103 }
104
105 G4DataSet::~G4DataSet()
106 {
107 delete algorithm;
108 if (energies) delete energies;
109 if (data) delete data;
110 if (pdf) delete pdf;
111 }
112
113 G4double G4DataSet::FindValue(G4double energy,
114 {
115 if (!energies) G4Exception("G4DataSet::FindV
116 "pii00000120",
117 FatalException,
118 "energies == 0");
119 if (energies->empty()) return 0;
120 if (energy <= (*energies)[0]) return (*data)
121
122 std::size_t i = energies->size()-1;
123 if (energy >= (*energies)[i]) return (*data)
124
125 G4double interpolated = algorithm->Calculate
126 return interpolated;
127 }
128
129
130 void G4DataSet::PrintData(void) const
131 {
132 if (!energies)
133 {
134 G4cout << "Data not available." << G4end
135 }
136 else
137 {
138 std::size_t size = energies->size();
139 for (std::size_t i(0); i<size; i++)
140 {
141 G4cout << "Point: " << ((*energies)[i]/uni
142 << " - Data value: " << ((*data)[i]/unitD
143 if (pdf != 0) G4cout << " - PDF : " << (*p
144 G4cout << G4endl;
145 }
146 }
147 }
148
149
150 void G4DataSet::SetEnergiesData(G4DataVector*
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",
163 FatalException,
164 "different size for energies and data (zer
165
166 if (energies == 0) return;
167
168 if (energies->size() != data->size())
169 G4Exception("G4DataSet::SetEnergiesData",
170 "pii00000131",
171 FatalException,
172 "different size for energies and data");
173 }
174
175 G4bool G4DataSet::LoadData(const G4String& fil
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
181 // -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 f
191
192 G4Exception("G4CompositeDataSet::LoadDat
193 "pii00000140",
194 FatalException,
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 =
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
228 G4bool G4DataSet::SaveData(const G4String& nam
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
234 // -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 - canno
244
245 G4Exception("G4CompositeDataSet::SaveDat
246 "pii00000150",
247 FatalException,
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-
259 G4DataVector::const_iterator endI(energi
260 G4DataVector::const_iterator j(data->beg
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
302 std::size_t G4DataSet::FindLowerBound(G4double
303 {
304 std::size_t lowerBound = 0;
305 std::size_t upperBound(energies->size() - 1)
306
307 while (lowerBound <= upperBound)
308 {
309 std::size_t midBin((lowerBound + upperBo
310
311 if (x < (*energies)[midBin]) upperBound
312 else lowerBound = midBin + 1;
313 }
314
315 return upperBound;
316 }
317
318
319 std::size_t G4DataSet::FindLowerBound(G4double
320 {
321 std::size_t lowerBound = 0;;
322 std::size_t upperBound(values->size() - 1);
323
324 while (lowerBound <= upperBound)
325 {
326 std::size_t midBin((lowerBound + upperBo
327
328 if (x < (*values)[midBin]) upperBound =
329 else lowerBound = midBin + 1;
330 }
331
332 return upperBound;
333 }
334
335
336 G4String G4DataSet::FullFileName(const G4Strin
337 {
338 const char* path = G4FindDataDir("G4PIIDATA"
339 if (!path)
340 G4Exception("G4DataSet::FullFileName",
341 "pii00000160",
342 FatalException,
343 "G4PIIDATA environment variable not set");
344
345 std::ostringstream fullFileName;
346 fullFileName << path << '/' << name << z <<
347
348 return G4String(fullFileName.str().c_str());
349 }
350
351
352 void G4DataSet::BuildPdf()
353 {
354 pdf = new G4DataVector;
355 G4Integrator <G4DataSet, G4double(G4DataSet:
356
357 std::size_t nData = data->size();
358 pdf->push_back(0.);
359
360 // Integrate the data distribution
361 std::size_t 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(thi
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 /* comp
383 {
384 // Random select a X value according to the
385 // derived from the data
386
387 if (!pdf) G4Exception("G4DataSet::RandomSele
388 "pii00000170",
389 FatalException,
390 "PDF has not been created for this data
391
392 G4double value = 0.;
393 G4double x = G4UniformRand();
394
395 // Locate the random value in the X vector b
396 G4int bin = (G4int)FindLowerBound(x,pdf);
397
398 // Interpolate the PDF to calculate the X va
399 // linear interpolation in the first bin (to
400 // interpolation with associated data set al
401
402 G4LinInterpolation linearAlgo;
403 if (bin == 0) value = linearAlgo.Calculate(x
404 else value = algorithm->Calculate(x, bin, *p
405
406 // G4cout << x << " random bin "<< bin << "
407 return value;
408 }
409
410 G4double G4DataSet::IntegrationFunction(G4doub
411 {
412 // This function is needed by G4Integrator t
413
414 G4double y = 0;
415
416 // Locate the random value in the X vector b
417 G4int bin = (G4int)FindLowerBound(x);
418
419 // Interpolate to calculate the X value:
420 // linear interpolation in the first bin (to
421 // interpolation with associated algorithm i
422
423 G4LinInterpolation linearAlgo;
424
425 if (bin == 0) y = linearAlgo.Calculate(x, bi
426 else y = algorithm->Calculate(x, bin, *energ
427
428 return y;
429 }
430
431