| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id: G4DataSet.cc 70904 2013-06-07 10:34:25Z gcosmo $
27 // 28 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 29 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 30 //
30 // History: 31 // History:
31 // ----------- 32 // -----------
32 // 31 Jul 2001 MGP Created 33 // 31 Jul 2001 MGP Created
33 // 31 Jul 2008 MGP Revised and rename 34 // 31 Jul 2008 MGP Revised and renamed to G4DataSet
34 // 35 //
35 // ------------------------------------------- 36 // -------------------------------------------------------------------
36 37
37 #include "G4DataSet.hh" 38 #include "G4DataSet.hh"
38 #include "G4IInterpolator.hh" 39 #include "G4IInterpolator.hh"
39 #include <fstream> 40 #include <fstream>
40 #include <sstream> 41 #include <sstream>
41 #include "G4Integrator.hh" 42 #include "G4Integrator.hh"
42 #include "Randomize.hh" 43 #include "Randomize.hh"
43 #include "G4LinInterpolation.hh" 44 #include "G4LinInterpolation.hh"
44 45
45 46
46 G4DataSet::G4DataSet(G4int Z, 47 G4DataSet::G4DataSet(G4int Z,
47 G4IInterpolator* algo, 48 G4IInterpolator* algo,
48 G4double xUnit, 49 G4double xUnit,
49 G4double yUnit, 50 G4double yUnit,
50 G4bool random): 51 G4bool random):
51 z(Z), 52 z(Z),
52 energies(0), 53 energies(0),
53 data(0), 54 data(0),
54 algorithm(algo), 55 algorithm(algo),
55 unitEnergies(xUnit), 56 unitEnergies(xUnit),
56 unitData(yUnit), 57 unitData(yUnit),
57 pdf(0), 58 pdf(0),
58 randomSet(random) 59 randomSet(random)
59 { 60 {
60 if (algorithm == 0) G4Exception("G4DataSet:: 61 if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
61 "pii00000101", 62 "pii00000101",
62 FatalException, 63 FatalException,
63 "Interpolation == 0"); 64 "Interpolation == 0");
64 if (randomSet) BuildPdf(); 65 if (randomSet) BuildPdf();
65 } 66 }
66 67
67 G4DataSet::G4DataSet(G4int argZ, 68 G4DataSet::G4DataSet(G4int argZ,
68 G4DataVector* dataX, 69 G4DataVector* dataX,
69 G4DataVector* dataY, 70 G4DataVector* dataY,
70 G4IInterpolator* algo, 71 G4IInterpolator* algo,
71 G4double xUnit, 72 G4double xUnit,
72 G4double yUnit, 73 G4double yUnit,
73 G4bool random): 74 G4bool random):
74 z(argZ), 75 z(argZ),
75 energies(dataX), 76 energies(dataX),
76 data(dataY), 77 data(dataY),
77 algorithm(algo), 78 algorithm(algo),
78 unitEnergies(xUnit), 79 unitEnergies(xUnit),
79 unitData(yUnit), 80 unitData(yUnit),
80 pdf(0), 81 pdf(0),
81 randomSet(random) 82 randomSet(random)
82 { 83 {
83 if (algorithm == 0) G4Exception("G4DataSet:: 84 if (algorithm == 0) G4Exception("G4DataSet::G4DataSet",
84 "pii00000110", 85 "pii00000110",
85 FatalException, 86 FatalException,
86 "Interpolation == 0"); 87 "Interpolation == 0");
87 88
88 if ((energies == 0) ^ (data == 0)) 89 if ((energies == 0) ^ (data == 0))
89 G4Exception("G4DataSet::G4DataSet", 90 G4Exception("G4DataSet::G4DataSet",
90 "pii00000111-", 91 "pii00000111-",
91 FatalException, 92 FatalException,
92 "different size for energies and data (zer 93 "different size for energies and data (zero case)");
93 94
94 if (energies == 0) return; 95 if (energies == 0) return;
95 96
96 if (energies->size() != data->size()) 97 if (energies->size() != data->size())
97 G4Exception("G4DataSet::G4DataSet", 98 G4Exception("G4DataSet::G4DataSet",
98 "pii00000112", 99 "pii00000112",
99 FatalException, 100 FatalException,
100 "different size for energies and data"); 101 "different size for energies and data");
101 102
102 if (randomSet) BuildPdf(); 103 if (randomSet) BuildPdf();
103 } 104 }
104 105
105 G4DataSet::~G4DataSet() 106 G4DataSet::~G4DataSet()
106 { 107 {
107 delete algorithm; 108 delete algorithm;
108 if (energies) delete energies; 109 if (energies) delete energies;
109 if (data) delete data; 110 if (data) delete data;
110 if (pdf) delete pdf; 111 if (pdf) delete pdf;
111 } 112 }
112 113
113 G4double G4DataSet::FindValue(G4double energy, 114 G4double G4DataSet::FindValue(G4double energy, G4int /* componentId */) const
114 { 115 {
115 if (!energies) G4Exception("G4DataSet::FindV 116 if (!energies) G4Exception("G4DataSet::FindValue",
116 "pii00000120", 117 "pii00000120",
117 FatalException, 118 FatalException,
118 "energies == 0"); 119 "energies == 0");
119 if (energies->empty()) return 0; 120 if (energies->empty()) return 0;
120 if (energy <= (*energies)[0]) return (*data) 121 if (energy <= (*energies)[0]) return (*data)[0];
121 122
122 std::size_t i = energies->size()-1; << 123 size_t i = energies->size()-1;
123 if (energy >= (*energies)[i]) return (*data) 124 if (energy >= (*energies)[i]) return (*data)[i];
124 125
125 G4double interpolated = algorithm->Calculate << 126 G4double interpolated = algorithm->Calculate(energy,FindLowerBound(energy),*energies,*data);
126 return interpolated; 127 return interpolated;
127 } 128 }
128 129
129 130
130 void G4DataSet::PrintData(void) const 131 void G4DataSet::PrintData(void) const
131 { 132 {
132 if (!energies) 133 if (!energies)
133 { 134 {
134 G4cout << "Data not available." << G4end 135 G4cout << "Data not available." << G4endl;
135 } 136 }
136 else 137 else
137 { 138 {
138 std::size_t size = energies->size(); << 139 size_t size = energies->size();
139 for (std::size_t i(0); i<size; i++) << 140 for (size_t i(0); i<size; i++)
140 { 141 {
141 G4cout << "Point: " << ((*energies)[i]/uni 142 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
142 << " - Data value: " << ((*data)[i]/unitD 143 << " - Data value: " << ((*data)[i]/unitData);
143 if (pdf != 0) G4cout << " - PDF : " << (*p 144 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
144 G4cout << G4endl; 145 G4cout << G4endl;
145 } 146 }
146 } 147 }
147 } 148 }
148 149
149 150
150 void G4DataSet::SetEnergiesData(G4DataVector* 151 void G4DataSet::SetEnergiesData(G4DataVector* dataX,
151 G4DataVector* dataY, 152 G4DataVector* dataY,
152 G4int /* componentId */) 153 G4int /* componentId */)
153 { 154 {
154 if (energies) delete energies; 155 if (energies) delete energies;
155 energies = dataX; 156 energies = dataX;
156 157
157 if (data) delete data; 158 if (data) delete data;
158 data = dataY; 159 data = dataY;
159 160
160 if ((energies == 0) ^ (data==0)) 161 if ((energies == 0) ^ (data==0))
161 G4Exception("G4DataSet::SetEnergiesData", 162 G4Exception("G4DataSet::SetEnergiesData",
162 "pii00000130", 163 "pii00000130",
163 FatalException, 164 FatalException,
164 "different size for energies and data (zer 165 "different size for energies and data (zero case)");
165 166
166 if (energies == 0) return; 167 if (energies == 0) return;
167 168
168 if (energies->size() != data->size()) 169 if (energies->size() != data->size())
169 G4Exception("G4DataSet::SetEnergiesData", 170 G4Exception("G4DataSet::SetEnergiesData",
170 "pii00000131", 171 "pii00000131",
171 FatalException, 172 FatalException,
172 "different size for energies and data"); 173 "different size for energies and data");
173 } 174 }
174 175
175 G4bool G4DataSet::LoadData(const G4String& fil 176 G4bool G4DataSet::LoadData(const G4String& fileName)
176 { 177 {
177 // The file is organized into two columns: 178 // The file is organized into two columns:
178 // 1st column is the energy 179 // 1st column is the energy
179 // 2nd column is the corresponding value 180 // 2nd column is the corresponding value
180 // The file terminates with the pattern: -1 181 // The file terminates with the pattern: -1 -1
181 // -2 182 // -2 -2
182 183
183 G4String fullFileName(FullFileName(fileName) 184 G4String fullFileName(FullFileName(fileName));
184 std::ifstream in(fullFileName); 185 std::ifstream in(fullFileName);
185 186
186 if (!in.is_open()) 187 if (!in.is_open())
187 { 188 {
188 189
189 std::ostringstream message; 190 std::ostringstream message;
190 message << "G4DataSet::LoadData - data f 191 message << "G4DataSet::LoadData - data file " << fullFileName << " not found";
191 192
192 G4Exception("G4CompositeDataSet::LoadDat 193 G4Exception("G4CompositeDataSet::LoadData",
193 "pii00000140", 194 "pii00000140",
194 FatalException, 195 FatalException,
195 message.str().c_str()); 196 message.str().c_str());
196 } 197 }
197 198
198 G4DataVector* argEnergies=new G4DataVector; 199 G4DataVector* argEnergies=new G4DataVector;
199 G4DataVector* argData=new G4DataVector; 200 G4DataVector* argData=new G4DataVector;
200 201
201 G4double a; 202 G4double a;
202 bool energyColumn(true); 203 bool energyColumn(true);
203 204
204 do 205 do
205 { 206 {
206 in >> a; 207 in >> a;
207 208
208 if (a!=-1 && a!=-2) 209 if (a!=-1 && a!=-2)
209 { 210 {
210 if (energyColumn) 211 if (energyColumn)
211 { 212 {
212 // std::cout << fullFileName << ", a = 213 // std::cout << fullFileName << ", a = " << a <<std::endl;
213 argEnergies->push_back(a*unitEnergies) 214 argEnergies->push_back(a*unitEnergies);
214 } 215 }
215 else 216 else
216 argData->push_back(a*unitData); 217 argData->push_back(a*unitData);
217 energyColumn=(!energyColumn); 218 energyColumn=(!energyColumn);
218 } 219 }
219 } 220 }
220 while (a != -2); 221 while (a != -2);
221 222
222 SetEnergiesData(argEnergies, argData, 0); 223 SetEnergiesData(argEnergies, argData, 0);
223 if (randomSet) BuildPdf(); 224 if (randomSet) BuildPdf();
224 225
225 return true; 226 return true;
226 } 227 }
227 228
228 G4bool G4DataSet::SaveData(const G4String& nam 229 G4bool G4DataSet::SaveData(const G4String& name) const
229 { 230 {
230 // The file is organized into two columns: 231 // The file is organized into two columns:
231 // 1st column is the energy 232 // 1st column is the energy
232 // 2nd column is the corresponding value 233 // 2nd column is the corresponding value
233 // The file terminates with the pattern: -1 234 // The file terminates with the pattern: -1 -1
234 // -2 235 // -2 -2
235 236
236 G4String fullFileName(FullFileName(name)); 237 G4String fullFileName(FullFileName(name));
237 std::ofstream out(fullFileName); 238 std::ofstream out(fullFileName);
238 239
239 if (!out.is_open()) 240 if (!out.is_open())
240 { 241 {
241 242
242 std::ostringstream message; 243 std::ostringstream message;
243 message << "G4DataSet:: SaveData - canno 244 message << "G4DataSet:: SaveData - cannot open " << fullFileName;
244 245
245 G4Exception("G4CompositeDataSet::SaveDat 246 G4Exception("G4CompositeDataSet::SaveData",
246 "pii00000150", 247 "pii00000150",
247 FatalException, 248 FatalException,
248 message.str().c_str()); 249 message.str().c_str());
249 250
250 } 251 }
251 252
252 out.precision(10); 253 out.precision(10);
253 out.width(15); 254 out.width(15);
254 out.setf(std::ofstream::left); 255 out.setf(std::ofstream::left);
255 256
256 if (energies!=0 && data!=0) 257 if (energies!=0 && data!=0)
257 { 258 {
258 G4DataVector::const_iterator i(energies- 259 G4DataVector::const_iterator i(energies->begin());
259 G4DataVector::const_iterator endI(energi 260 G4DataVector::const_iterator endI(energies->end());
260 G4DataVector::const_iterator j(data->beg 261 G4DataVector::const_iterator j(data->begin());
261 262
262 while (i!=endI) 263 while (i!=endI)
263 { 264 {
264 out.precision(10); 265 out.precision(10);
265 out.width(15); 266 out.width(15);
266 out.setf(std::ofstream::left); 267 out.setf(std::ofstream::left);
267 out << ((*i)/unitEnergies) << ' '; 268 out << ((*i)/unitEnergies) << ' ';
268 269
269 out.precision(10); 270 out.precision(10);
270 out.width(15); 271 out.width(15);
271 out.setf(std::ofstream::left); 272 out.setf(std::ofstream::left);
272 out << ((*j)/unitData) << std::endl; 273 out << ((*j)/unitData) << std::endl;
273 274
274 i++; 275 i++;
275 j++; 276 j++;
276 } 277 }
277 } 278 }
278 279
279 out.precision(10); 280 out.precision(10);
280 out.width(15); 281 out.width(15);
281 out.setf(std::ofstream::left); 282 out.setf(std::ofstream::left);
282 out << -1.f << ' '; 283 out << -1.f << ' ';
283 284
284 out.precision(10); 285 out.precision(10);
285 out.width(15); 286 out.width(15);
286 out.setf(std::ofstream::left); 287 out.setf(std::ofstream::left);
287 out << -1.f << std::endl; 288 out << -1.f << std::endl;
288 289
289 out.precision(10); 290 out.precision(10);
290 out.width(15); 291 out.width(15);
291 out.setf(std::ofstream::left); 292 out.setf(std::ofstream::left);
292 out << -2.f << ' '; 293 out << -2.f << ' ';
293 294
294 out.precision(10); 295 out.precision(10);
295 out.width(15); 296 out.width(15);
296 out.setf(std::ofstream::left); 297 out.setf(std::ofstream::left);
297 out << -2.f << std::endl; 298 out << -2.f << std::endl;
298 299
299 return true; 300 return true;
300 } 301 }
301 302
302 std::size_t G4DataSet::FindLowerBound(G4double << 303 size_t G4DataSet::FindLowerBound(G4double x) const
303 { 304 {
304 std::size_t lowerBound = 0; << 305 size_t lowerBound = 0;
305 std::size_t upperBound(energies->size() - 1) << 306 size_t upperBound(energies->size() - 1);
306 307
307 while (lowerBound <= upperBound) 308 while (lowerBound <= upperBound)
308 { 309 {
309 std::size_t midBin((lowerBound + upperBo << 310 size_t midBin((lowerBound + upperBound) / 2);
310 311
311 if (x < (*energies)[midBin]) upperBound 312 if (x < (*energies)[midBin]) upperBound = midBin - 1;
312 else lowerBound = midBin + 1; 313 else lowerBound = midBin + 1;
313 } 314 }
314 315
315 return upperBound; 316 return upperBound;
316 } 317 }
317 318
318 319
319 std::size_t G4DataSet::FindLowerBound(G4double << 320 size_t G4DataSet::FindLowerBound(G4double x, G4DataVector* values) const
320 { 321 {
321 std::size_t lowerBound = 0;; << 322 size_t lowerBound = 0;;
322 std::size_t upperBound(values->size() - 1); << 323 size_t upperBound(values->size() - 1);
323 324
324 while (lowerBound <= upperBound) 325 while (lowerBound <= upperBound)
325 { 326 {
326 std::size_t midBin((lowerBound + upperBo << 327 size_t midBin((lowerBound + upperBound) / 2);
327 328
328 if (x < (*values)[midBin]) upperBound = 329 if (x < (*values)[midBin]) upperBound = midBin - 1;
329 else lowerBound = midBin + 1; 330 else lowerBound = midBin + 1;
330 } 331 }
331 332
332 return upperBound; 333 return upperBound;
333 } 334 }
334 335
335 336
336 G4String G4DataSet::FullFileName(const G4Strin 337 G4String G4DataSet::FullFileName(const G4String& name) const
337 { 338 {
338 const char* path = G4FindDataDir("G4PIIDATA" << 339 char* path = getenv("G4PIIDATA");
339 if (!path) 340 if (!path)
340 G4Exception("G4DataSet::FullFileName", 341 G4Exception("G4DataSet::FullFileName",
341 "pii00000160", 342 "pii00000160",
342 FatalException, 343 FatalException,
343 "G4PIIDATA environment variable not set"); 344 "G4PIIDATA environment variable not set");
344 345
345 std::ostringstream fullFileName; 346 std::ostringstream fullFileName;
346 fullFileName << path << '/' << name << z << 347 fullFileName << path << '/' << name << z << ".dat";
347 348
348 return G4String(fullFileName.str().c_str()); 349 return G4String(fullFileName.str().c_str());
349 } 350 }
350 351
351 352
352 void G4DataSet::BuildPdf() 353 void G4DataSet::BuildPdf()
353 { 354 {
354 pdf = new G4DataVector; 355 pdf = new G4DataVector;
355 G4Integrator <G4DataSet, G4double(G4DataSet: 356 G4Integrator <G4DataSet, G4double(G4DataSet::*)(G4double)> integrator;
356 357
357 std::size_t nData = data->size(); << 358 G4int nData = data->size();
358 pdf->push_back(0.); 359 pdf->push_back(0.);
359 360
360 // Integrate the data distribution 361 // Integrate the data distribution
361 std::size_t i; << 362 G4int i;
362 G4double totalSum = 0.; 363 G4double totalSum = 0.;
363 for (i=1; i<nData; i++) 364 for (i=1; i<nData; i++)
364 { 365 {
365 G4double xLow = (*energies)[i-1]; 366 G4double xLow = (*energies)[i-1];
366 G4double xHigh = (*energies)[i]; 367 G4double xHigh = (*energies)[i];
367 G4double sum = integrator.Legendre96(thi 368 G4double sum = integrator.Legendre96(this, &G4DataSet::IntegrationFunction, xLow, xHigh);
368 totalSum = totalSum + sum; 369 totalSum = totalSum + sum;
369 pdf->push_back(totalSum); 370 pdf->push_back(totalSum);
370 } 371 }
371 372
372 // Normalize to the last bin 373 // Normalize to the last bin
373 G4double tot = 0.; 374 G4double tot = 0.;
374 if (totalSum > 0.) tot = 1. / totalSum; 375 if (totalSum > 0.) tot = 1. / totalSum;
375 for (i=1; i<nData; i++) 376 for (i=1; i<nData; i++)
376 { 377 {
377 (*pdf)[i] = (*pdf)[i] * tot; 378 (*pdf)[i] = (*pdf)[i] * tot;
378 } 379 }
379 } 380 }
380 381
381 382
382 G4double G4DataSet::RandomSelect(G4int /* comp 383 G4double G4DataSet::RandomSelect(G4int /* componentId */) const
383 { 384 {
384 // Random select a X value according to the 385 // Random select a X value according to the cumulative probability distribution
385 // derived from the data 386 // derived from the data
386 387
387 if (!pdf) G4Exception("G4DataSet::RandomSele 388 if (!pdf) G4Exception("G4DataSet::RandomSelect",
388 "pii00000170", 389 "pii00000170",
389 FatalException, 390 FatalException,
390 "PDF has not been created for this data 391 "PDF has not been created for this data set");
391 392
392 G4double value = 0.; 393 G4double value = 0.;
393 G4double x = G4UniformRand(); 394 G4double x = G4UniformRand();
394 395
395 // Locate the random value in the X vector b 396 // Locate the random value in the X vector based on the PDF
396 G4int bin = (G4int)FindLowerBound(x,pdf); << 397 size_t bin = FindLowerBound(x,pdf);
397 398
398 // Interpolate the PDF to calculate the X va 399 // Interpolate the PDF to calculate the X value:
399 // linear interpolation in the first bin (to 400 // linear interpolation in the first bin (to avoid problem with 0),
400 // interpolation with associated data set al 401 // interpolation with associated data set algorithm in other bins
401 402
402 G4LinInterpolation linearAlgo; 403 G4LinInterpolation linearAlgo;
403 if (bin == 0) value = linearAlgo.Calculate(x 404 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
404 else value = algorithm->Calculate(x, bin, *p 405 else value = algorithm->Calculate(x, bin, *pdf, *energies);
405 406
406 // G4cout << x << " random bin "<< bin << " 407 // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
407 return value; 408 return value;
408 } 409 }
409 410
410 G4double G4DataSet::IntegrationFunction(G4doub 411 G4double G4DataSet::IntegrationFunction(G4double x)
411 { 412 {
412 // This function is needed by G4Integrator t 413 // This function is needed by G4Integrator to calculate the integral of the data distribution
413 414
414 G4double y = 0; 415 G4double y = 0;
415 416
416 // Locate the random value in the X vector b 417 // Locate the random value in the X vector based on the PDF
417 G4int bin = (G4int)FindLowerBound(x); << 418 size_t bin = FindLowerBound(x);
418 419
419 // Interpolate to calculate the X value: 420 // Interpolate to calculate the X value:
420 // linear interpolation in the first bin (to 421 // linear interpolation in the first bin (to avoid problem with 0),
421 // interpolation with associated algorithm i 422 // interpolation with associated algorithm in other bins
422 423
423 G4LinInterpolation linearAlgo; 424 G4LinInterpolation linearAlgo;
424 425
425 if (bin == 0) y = linearAlgo.Calculate(x, bi 426 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
426 else y = algorithm->Calculate(x, bin, *energ 427 else y = algorithm->Calculate(x, bin, *energies, *data);
427 428
428 return y; 429 return y;
429 } 430 }
430 431
431 432