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