source: trunk/source/processes/electromagnetic/pii/src/G4DataSet.cc@ 1357

Last change on this file since 1357 was 1350, checked in by garnier, 15 years ago

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