| 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 // History: <<
27 // ----------- <<
28 // 31 Jul 2001 MGP Created <<
29 // <<
30 // 15 Jul 2009 Nicolas A. Karakatsanis <<
31 // <<
32 // - New Constructor <<
33 // to enhance CPU <<
34 // <<
35 // - Destructor was <<
36 // 26 //
37 // - LoadNonLogData << 27 // $Id: G4EMDataSet.cc,v 1.18 2008/03/17 13:40:53 pia Exp $
38 // dataset. It is << 28 // GEANT4 tag $Name: geant4-09-02 $
39 // 29 //
40 // - LoadData method << 30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
41 // It retrieves th <<
42 // respective log <<
43 // 31 //
44 // - SetLogEnergiesD << 32 // History:
45 // The EM data set << 33 // -----------
46 // These initializ << 34 // 31 Jul 2001 MGP Created
47 // operations <<
48 // <<
49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war <<
50 // 35 //
51 // ------------------------------------------- 36 // -------------------------------------------------------------------
52 37
53 #include "G4EMDataSet.hh" 38 #include "G4EMDataSet.hh"
54 #include "G4VDataSetAlgorithm.hh" 39 #include "G4VDataSetAlgorithm.hh"
55 #include <fstream> 40 #include <fstream>
56 #include <sstream> 41 #include <sstream>
57 #include "G4Integrator.hh" 42 #include "G4Integrator.hh"
58 #include "Randomize.hh" 43 #include "Randomize.hh"
59 #include "G4LinInterpolation.hh" 44 #include "G4LinInterpolation.hh"
60 45
>> 46
61 G4EMDataSet::G4EMDataSet(G4int Z, 47 G4EMDataSet::G4EMDataSet(G4int Z,
62 G4VDataSetAlgorithm* algo, 48 G4VDataSetAlgorithm* algo,
63 G4double xUnit, 49 G4double xUnit,
64 G4double yUnit, 50 G4double yUnit,
65 G4bool random): << 51 G4bool random):
66 energies(nullptr), << 52 z(Z),
67 data(nullptr), << 53 energies(0),
68 log_energies(nullptr), << 54 data(0),
69 log_data(nullptr), <<
70 algorithm(algo), 55 algorithm(algo),
71 pdf(nullptr), <<
72 unitEnergies(xUnit), 56 unitEnergies(xUnit),
73 unitData(yUnit), << 57 unitData(yUnit),
74 z(Z), << 58 pdf(0),
75 randomSet(random) 59 randomSet(random)
76 { 60 {
77 if (algorithm == nullptr) { << 61 if (algorithm == 0) G4Exception("G4EMDataSet::G4EMDataSet - interpolation == 0");
78 G4Exception("G4EMDataSet::G4EMDataSet", << 62 if (randomSet) BuildPdf();
79 "em1012",FatalException,"interpolation <<
80 } else if (randomSet) { BuildPdf(); } <<
81 } 63 }
82 64
83 //....oooOO0OOooo........oooOO0OOooo........oo <<
84 <<
85 G4EMDataSet::G4EMDataSet(G4int argZ, 65 G4EMDataSet::G4EMDataSet(G4int argZ,
86 G4DataVector* dataX, 66 G4DataVector* dataX,
87 G4DataVector* dataY, 67 G4DataVector* dataY,
88 G4VDataSetAlgorithm* algo, 68 G4VDataSetAlgorithm* algo,
89 G4double xUnit, 69 G4double xUnit,
90 G4double yUnit, 70 G4double yUnit,
91 G4bool random): 71 G4bool random):
92 energies(dataX), <<
93 data(dataY), <<
94 log_energies(nullptr), <<
95 log_data(nullptr), <<
96 algorithm(algo), <<
97 pdf(nullptr), <<
98 unitEnergies(xUnit), <<
99 unitData(yUnit), <<
100 z(argZ), 72 z(argZ),
101 randomSet(random) <<
102 { <<
103 if (!algorithm || !energies || !data) { <<
104 G4Exception("G4EMDataSet::G4EMDataSet", <<
105 "em1012",FatalException,"interpolation <<
106 } else { <<
107 if (energies->size() != data->size()) { <<
108 G4Exception("G4EMDataSet::G4EMDataSet", <<
109 "em1012",FatalException,"different size <<
110 } else if (randomSet) { <<
111 BuildPdf(); <<
112 } <<
113 } <<
114 } <<
115 <<
116 //....oooOO0OOooo........oooOO0OOooo........oo <<
117 <<
118 G4EMDataSet::G4EMDataSet(G4int argZ, <<
119 G4DataVector* dataX, <<
120 G4DataVector* dataY, <<
121 G4DataVector* dataLog <<
122 G4DataVector* dataLog <<
123 G4VDataSetAlgorithm* algo, <<
124 G4double xUnit, <<
125 G4double yUnit, <<
126 G4bool random): <<
127 energies(dataX), 73 energies(dataX),
128 data(dataY), 74 data(dataY),
129 log_energies(dataLogX), <<
130 log_data(dataLogY), <<
131 algorithm(algo), 75 algorithm(algo),
132 pdf(nullptr), <<
133 unitEnergies(xUnit), 76 unitEnergies(xUnit),
134 unitData(yUnit), 77 unitData(yUnit),
135 z(argZ), << 78 pdf(0),
136 randomSet(random) 79 randomSet(random)
137 { 80 {
138 if (!algorithm || !energies || !data || !log << 81 if (algorithm == 0) G4Exception("G4EMDataSet::G4EMDataSet - interpolation == 0");
139 G4Exception("G4EMDataSet::G4EMDataSet", <<
140 "em1012",FatalException,"interpolation <<
141 } else { <<
142 if ((energies->size() != data->size()) || <<
143 (energies->size() != log_energies->size()) | <<
144 (energies->size() != log_data->size())) { <<
145 G4Exception("G4EMDataSet::G4EMDataSet", <<
146 "em1012",FatalException,"different size <<
147 } else if (randomSet) { <<
148 BuildPdf(); <<
149 } <<
150 } <<
151 } <<
152 82
153 //....oooOO0OOooo........oooOO0OOooo........oo << 83 if ((energies == 0) ^ (data == 0))
>> 84 G4Exception("G4EMDataSet::G4EMDataSet - different size for energies and data (zero case)");
>> 85
>> 86 if (energies == 0) return;
>> 87
>> 88 if (energies->size() != data->size())
>> 89 G4Exception("G4EMDataSet::G4EMDataSet - different size for energies and data");
>> 90
>> 91 if (randomSet) BuildPdf();
>> 92 }
154 93
155 G4EMDataSet::~G4EMDataSet() 94 G4EMDataSet::~G4EMDataSet()
156 { 95 {
157 delete algorithm; 96 delete algorithm;
158 delete energies; << 97 if (energies) delete energies;
159 delete data; << 98 if (data) delete data;
160 delete pdf; << 99 if (pdf) delete pdf;
161 delete log_energies; <<
162 delete log_data; <<
163 } 100 }
164 101
165 //....oooOO0OOooo........oooOO0OOooo........oo <<
166 <<
167 G4double G4EMDataSet::FindValue(G4double energ 102 G4double G4EMDataSet::FindValue(G4double energy, G4int /* componentId */) const
168 { 103 {
169 if (energy <= (*energies)[0]) { << 104 if (!energies) G4Exception("G4EMDataSet::FindValue - energies == 0");
170 return (*data)[0]; << 105 if (energies->empty()) return 0;
171 } << 106 if (energy <= (*energies)[0]) return (*data)[0];
172 std::size_t i = energies->size()-1; << 107
173 if (energy >= (*energies)[i]) { return (*dat << 108 size_t i = energies->size()-1;
174 << 109 if (energy >= (*energies)[i]) return (*data)[i];
175 //Nicolas A. Karakatsanis: Check if the loga <<
176 // which Interpolati <<
177 if (log_energies != nullptr) <<
178 { <<
179 return algorithm->Calculate(energy, (G4in <<
180 } <<
181 110
182 return algorithm->Calculate(energy, (G4int)F << 111 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
183 } 112 }
184 113
185 //....oooOO0OOooo........oooOO0OOooo........oo <<
186 114
187 void G4EMDataSet::PrintData(void) const 115 void G4EMDataSet::PrintData(void) const
188 { 116 {
189 std::size_t size = energies->size(); << 117 if (!energies)
190 for (std::size_t i(0); i<size; ++i) <<
191 { 118 {
192 G4cout << "Point: " << ((*energies)[i]/u << 119 G4cout << "Data not available." << G4endl;
193 << " - Data value: " << ((*data)[i]/uni << 120 }
194 if (pdf != 0) G4cout << " - PDF : " << ( << 121 else
195 G4cout << G4endl; << 122 {
>> 123 size_t size = energies->size();
>> 124 for (size_t i(0); i<size; i++)
>> 125 {
>> 126 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
>> 127 << " - Data value: " << ((*data)[i]/unitData);
>> 128 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
>> 129 G4cout << G4endl;
>> 130 }
196 } 131 }
197 } 132 }
198 133
199 //....oooOO0OOooo........oooOO0OOooo........oo <<
200 134
201 void G4EMDataSet::SetEnergiesData(G4DataVector 135 void G4EMDataSet::SetEnergiesData(G4DataVector* dataX,
202 G4DataVector* dataY, 136 G4DataVector* dataY,
203 G4int /* componentId */) 137 G4int /* componentId */)
204 { 138 {
205 if(!dataX || !dataY) { << 139 if (energies) delete energies;
206 G4Exception("G4EMDataSet::SetEnergiesData" << 140 energies = dataX;
207 "em1012",FatalException,"new interpola <<
208 } else { <<
209 if (dataX->size() != dataY->size()) { <<
210 G4Exception("G4EMDataSet::SetEnergiesDat <<
211 "em1012",FatalException,"different size <<
212 } else { <<
213 <<
214 delete energies; <<
215 energies = dataX; <<
216 <<
217 delete data; <<
218 data = dataY; <<
219 //G4cout << "Size of energies: " << ener <<
220 //<< "Size of data: " << data->size() << <<
221 } <<
222 } <<
223 } <<
224 141
225 //....oooOO0OOooo........oooOO0OOooo........oo << 142 if (data) delete data;
>> 143 data = dataY;
>> 144
>> 145 if ((energies == 0) ^ (data==0))
>> 146 G4Exception("G4EMDataSet::SetEnergiesData - different size for energies and data (zero case)");
226 147
227 void G4EMDataSet::SetLogEnergiesData(G4DataVec << 148 if (energies == 0) return;
228 G4DataVec << 149
229 G4DataVec << 150 if (energies->size() != data->size())
230 G4DataVector* data_logY, << 151 G4Exception("G4EMDataSet::SetEnergiesData - different size for energies and data");
231 G4int /* <<
232 { <<
233 if(!dataX || !dataY || !data_logX || !data_l <<
234 G4Exception("G4EMDataSet::SetEnergiesData" <<
235 "em1012",FatalException,"new interpola <<
236 } else { <<
237 if (dataX->size() != dataY->size() || <<
238 dataX->size() != data_logX->size() || <<
239 dataX->size() != data_logY->size()) { <<
240 G4Exception("G4EMDataSet::SetEnergiesDat <<
241 "em1012",FatalException,"different size <<
242 } else { <<
243 <<
244 delete energies; <<
245 energies = dataX; <<
246 <<
247 delete data; <<
248 data = dataY; <<
249 <<
250 delete log_energies; <<
251 log_energies = data_logX; <<
252 <<
253 delete log_data; <<
254 log_data = data_logY; <<
255 //G4cout << "Size of energies: " << ener <<
256 //<< "Size of data: " << data->size() << <<
257 } <<
258 } <<
259 } 152 }
260 153
261 //....oooOO0OOooo........oooOO0OOooo........oo <<
262 <<
263 G4bool G4EMDataSet::LoadData(const G4String& f 154 G4bool G4EMDataSet::LoadData(const G4String& fileName)
264 { 155 {
265 // The file is organized into four columns: << 156 // The file is organized into two columns:
266 // 1st column contains the values of energy << 157 // 1st column is the energy
267 // 2nd column contains the corresponding data << 158 // 2nd column is the corresponding value
268 // The file terminates with the pattern: -1 << 159 // The file terminates with the pattern: -1 -1
269 // -2 << 160 // -2 -2
270 <<
271 G4String fullFileName(FullFileName(fileName) <<
272 std::ifstream in(fullFileName); <<
273 <<
274 if (!in.is_open()) <<
275 { <<
276 G4String message("data file \""); <<
277 message += fullFileName; <<
278 message += "\" not found"; <<
279 G4Exception("G4EMDataSet::LoadData", <<
280 "em1012",FatalException,message); <<
281 return false; <<
282 } <<
283 <<
284 delete energies; <<
285 delete data; <<
286 delete log_energies; <<
287 delete log_data; <<
288 energies = new G4DataVector; <<
289 data = new G4DataVector; <<
290 log_energies = new G4DataVector; <<
291 log_data = new G4DataVector; <<
292 <<
293 G4double a, b; <<
294 do <<
295 { <<
296 in >> a >> b; <<
297 <<
298 if (a != -1 && a != -2) <<
299 { <<
300 if (a==0.) { a=1e-300; } <<
301 if (b==0.) { b=1e-300; } <<
302 a *= unitEnergies; <<
303 b *= unitData; <<
304 energies->push_back(a); <<
305 log_energies->push_back(std::log10(a)); <<
306 data->push_back(b); <<
307 log_data->push_back(std::log10(b)); <<
308 } <<
309 } <<
310 while (a != -2); <<
311 <<
312 if (randomSet) { BuildPdf(); } <<
313 161
314 return true; <<
315 } <<
316 <<
317 //....oooOO0OOooo........oooOO0OOooo........oo <<
318 <<
319 G4bool G4EMDataSet::LoadNonLogData(const G4Str <<
320 { <<
321 // The file is organized into four columns: <<
322 // 1st column contains the values of energy <<
323 // 2nd column contains the corresponding data <<
324 // The file terminates with the pattern: -1 <<
325 // -2 <<
326 <<
327 G4String fullFileName(FullFileName(fileName) 162 G4String fullFileName(FullFileName(fileName));
328 std::ifstream in(fullFileName); 163 std::ifstream in(fullFileName);
>> 164
329 if (!in.is_open()) 165 if (!in.is_open())
330 { 166 {
331 G4String message("data file \""); << 167 G4String message("G4EMDataSet::LoadData - data file \"");
332 message += fullFileName; 168 message += fullFileName;
333 message += "\" not found"; 169 message += "\" not found";
334 G4Exception("G4EMDataSet::LoadNonLogData << 170 G4Exception(message);
335 "em1012",FatalException,message); <<
336 } 171 }
337 172
338 G4DataVector* argEnergies=new G4DataVector; 173 G4DataVector* argEnergies=new G4DataVector;
339 G4DataVector* argData=new G4DataVector; 174 G4DataVector* argData=new G4DataVector;
340 175
341 G4double a; 176 G4double a;
342 G4int k = 0; << 177 bool energyColumn(true);
343 G4int nColumns = 2; <<
344 178
345 do 179 do
346 { 180 {
347 in >> a; 181 in >> a;
348 182
349 if (a != -1 && a != -2) << 183 if (a!=-1 && a!=-2)
350 { 184 {
351 if (k%nColumns == 0) << 185 if (energyColumn)
352 { << 186 argEnergies->push_back(a*unitEnergies);
353 argEnergies->push_back(a*unitEnergies); << 187 else
354 } << 188 argData->push_back(a*unitData);
355 else if (k%nColumns == 1) << 189 energyColumn=(!energyColumn);
356 { <<
357 argData->push_back(a*unitData); <<
358 } <<
359 k++; <<
360 } 190 }
361 } 191 }
362 while (a != -2); 192 while (a != -2);
363 193
364 SetEnergiesData(argEnergies, argData, 0); 194 SetEnergiesData(argEnergies, argData, 0);
365 <<
366 if (randomSet) BuildPdf(); 195 if (randomSet) BuildPdf();
367 196
368 return true; 197 return true;
369 } 198 }
370 199
371 //....oooOO0OOooo........oooOO0OOooo........oo <<
372 <<
373 G4bool G4EMDataSet::SaveData(const G4String& n 200 G4bool G4EMDataSet::SaveData(const G4String& name) const
374 { 201 {
375 // The file is organized into two columns: 202 // The file is organized into two columns:
376 // 1st column is the energy 203 // 1st column is the energy
377 // 2nd column is the corresponding value 204 // 2nd column is the corresponding value
378 // The file terminates with the pattern: -1 205 // The file terminates with the pattern: -1 -1
379 // -2 206 // -2 -2
380 207
381 G4String fullFileName(FullFileName(name)); 208 G4String fullFileName(FullFileName(name));
382 std::ofstream out(fullFileName); 209 std::ofstream out(fullFileName);
383 210
384 if (!out.is_open()) 211 if (!out.is_open())
385 { 212 {
386 G4String message("cannot open \""); << 213 G4String message("G4EMDataSet::SaveData - cannot open \"");
387 message+=fullFileName; 214 message+=fullFileName;
388 message+="\""; 215 message+="\"";
389 G4Exception("G4EMDataSet::SaveData", << 216 G4Exception(message);
390 "em1012",FatalException,message); <<
391 } 217 }
392 218
393 out.precision(10); 219 out.precision(10);
394 out.width(15); 220 out.width(15);
395 out.setf(std::ofstream::left); 221 out.setf(std::ofstream::left);
396 222
397 if (energies!=0 && data!=0) 223 if (energies!=0 && data!=0)
398 { 224 {
399 G4DataVector::const_iterator i(energies- 225 G4DataVector::const_iterator i(energies->begin());
400 G4DataVector::const_iterator endI(energi 226 G4DataVector::const_iterator endI(energies->end());
401 G4DataVector::const_iterator j(data->beg 227 G4DataVector::const_iterator j(data->begin());
402 228
403 while (i!=endI) 229 while (i!=endI)
404 { 230 {
405 out.precision(10); 231 out.precision(10);
406 out.width(15); 232 out.width(15);
407 out.setf(std::ofstream::left); 233 out.setf(std::ofstream::left);
408 out << ((*i)/unitEnergies) << ' '; 234 out << ((*i)/unitEnergies) << ' ';
409 235
410 out.precision(10); 236 out.precision(10);
411 out.width(15); 237 out.width(15);
412 out.setf(std::ofstream::left); 238 out.setf(std::ofstream::left);
413 out << ((*j)/unitData) << std::endl; 239 out << ((*j)/unitData) << std::endl;
414 240
415 i++; 241 i++;
416 j++; 242 j++;
417 } 243 }
418 } 244 }
419 245
420 out.precision(10); 246 out.precision(10);
421 out.width(15); 247 out.width(15);
422 out.setf(std::ofstream::left); 248 out.setf(std::ofstream::left);
423 out << -1.f << ' '; 249 out << -1.f << ' ';
424 250
425 out.precision(10); 251 out.precision(10);
426 out.width(15); 252 out.width(15);
427 out.setf(std::ofstream::left); 253 out.setf(std::ofstream::left);
428 out << -1.f << std::endl; 254 out << -1.f << std::endl;
429 255
430 out.precision(10); 256 out.precision(10);
431 out.width(15); 257 out.width(15);
432 out.setf(std::ofstream::left); 258 out.setf(std::ofstream::left);
433 out << -2.f << ' '; 259 out << -2.f << ' ';
434 260
435 out.precision(10); 261 out.precision(10);
436 out.width(15); 262 out.width(15);
437 out.setf(std::ofstream::left); 263 out.setf(std::ofstream::left);
438 out << -2.f << std::endl; 264 out << -2.f << std::endl;
439 265
440 return true; 266 return true;
441 } 267 }
442 268
443 //....oooOO0OOooo........oooOO0OOooo........oo << 269 size_t G4EMDataSet::FindLowerBound(G4double x) const
444 <<
445 std::size_t G4EMDataSet::FindLowerBound(G4doub <<
446 { 270 {
447 std::size_t lowerBound = 0; << 271 size_t lowerBound = 0;
448 std::size_t upperBound(energies->size() - 1) << 272 size_t upperBound(energies->size() - 1);
449 273
450 while (lowerBound <= upperBound) 274 while (lowerBound <= upperBound)
451 { 275 {
452 std::size_t midBin((lowerBound + upperBo << 276 size_t midBin((lowerBound + upperBound) / 2);
453 277
454 if (x < (*energies)[midBin]) upperBound 278 if (x < (*energies)[midBin]) upperBound = midBin - 1;
455 else lowerBound = midBin + 1; 279 else lowerBound = midBin + 1;
456 } 280 }
457 281
458 return upperBound; 282 return upperBound;
459 } 283 }
460 284
461 //....oooOO0OOooo........oooOO0OOooo........oo <<
462 285
463 std::size_t G4EMDataSet::FindLowerBound(G4doub << 286 size_t G4EMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
464 { 287 {
465 std::size_t lowerBound = 0;; << 288 size_t lowerBound = 0;;
466 std::size_t upperBound(values->size() - 1); << 289 size_t upperBound(values->size() - 1);
467 290
468 while (lowerBound <= upperBound) 291 while (lowerBound <= upperBound)
469 { 292 {
470 std::size_t midBin((lowerBound + upperBo << 293 size_t midBin((lowerBound + upperBound) / 2);
471 294
472 if (x < (*values)[midBin]) upperBound = 295 if (x < (*values)[midBin]) upperBound = midBin - 1;
473 else lowerBound = midBin + 1; 296 else lowerBound = midBin + 1;
474 } 297 }
475 298
476 return upperBound; 299 return upperBound;
477 } 300 }
478 301
479 //....oooOO0OOooo........oooOO0OOooo........oo <<
480 302
481 G4String G4EMDataSet::FullFileName(const G4Str 303 G4String G4EMDataSet::FullFileName(const G4String& name) const
482 { 304 {
483 const char* path = G4FindDataDir("G4LEDATA") << 305 char* path = getenv("G4LEDATA");
484 if (!path) { << 306 if (!path)
485 G4Exception("G4EMDataSet::FullFileName", << 307 G4Exception("G4EMDataSet::FullFileName - G4LEDATA environment variable not set");
486 "em0006",FatalException,"G4LEDATA envi << 308
487 return ""; <<
488 } <<
489 std::ostringstream fullFileName; 309 std::ostringstream fullFileName;
490 fullFileName << path << '/' << name << z << 310 fullFileName << path << '/' << name << z << ".dat";
491 311
492 return G4String(fullFileName.str().c_str()); 312 return G4String(fullFileName.str().c_str());
493 } 313 }
494 314
495 //....oooOO0OOooo........oooOO0OOooo........oo <<
496 315
497 void G4EMDataSet::BuildPdf() 316 void G4EMDataSet::BuildPdf()
498 { 317 {
499 pdf = new G4DataVector; 318 pdf = new G4DataVector;
500 G4Integrator <G4EMDataSet, G4double(G4EMData 319 G4Integrator <G4EMDataSet, G4double(G4EMDataSet::*)(G4double)> integrator;
501 320
502 std::size_t nData = data->size(); << 321 G4int nData = data->size();
503 pdf->push_back(0.); 322 pdf->push_back(0.);
504 323
505 // Integrate the data distribution 324 // Integrate the data distribution
506 std::size_t i; << 325 G4int i;
507 G4double totalSum = 0.; 326 G4double totalSum = 0.;
508 for (i=1; i<nData; ++i) << 327 for (i=1; i<nData; i++)
509 { 328 {
510 G4double xLow = (*energies)[i-1]; 329 G4double xLow = (*energies)[i-1];
511 G4double xHigh = (*energies)[i]; 330 G4double xHigh = (*energies)[i];
512 G4double sum = integrator.Legendre96(thi 331 G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh);
513 totalSum = totalSum + sum; 332 totalSum = totalSum + sum;
514 pdf->push_back(totalSum); 333 pdf->push_back(totalSum);
515 } 334 }
516 335
517 // Normalize to the last bin 336 // Normalize to the last bin
518 G4double tot = 0.; 337 G4double tot = 0.;
519 if (totalSum > 0.) tot = 1. / totalSum; 338 if (totalSum > 0.) tot = 1. / totalSum;
520 for (i=1; i<nData; ++i) << 339 for (i=1; i<nData; i++)
521 { 340 {
522 (*pdf)[i] = (*pdf)[i] * tot; 341 (*pdf)[i] = (*pdf)[i] * tot;
523 } 342 }
524 } 343 }
525 344
526 //....oooOO0OOooo........oooOO0OOooo........oo <<
527 345
528 G4double G4EMDataSet::RandomSelect(G4int /* co 346 G4double G4EMDataSet::RandomSelect(G4int /* componentId */) const
529 { 347 {
530 G4double value = 0.; <<
531 // Random select a X value according to the 348 // Random select a X value according to the cumulative probability distribution
532 // derived from the data 349 // derived from the data
533 350
534 if (!pdf) { << 351 if (!pdf) G4Exception("G4EMDataSet::RandomSelect - PDF has not been created for this data set");
535 G4Exception("G4EMDataSet::RandomSelect", <<
536 "em1012",FatalException,"PDF has not b <<
537 return value; <<
538 } <<
539 352
>> 353 G4double value = 0.;
540 G4double x = G4UniformRand(); 354 G4double x = G4UniformRand();
541 355
542 // Locate the random value in the X vector b 356 // Locate the random value in the X vector based on the PDF
543 G4int bin = (G4int)FindLowerBound(x,pdf); << 357 size_t bin = FindLowerBound(x,pdf);
544 358
545 // Interpolate the PDF to calculate the X va 359 // Interpolate the PDF to calculate the X value:
546 // linear interpolation in the first bin (to 360 // linear interpolation in the first bin (to avoid problem with 0),
547 // interpolation with associated data set al 361 // interpolation with associated data set algorithm in other bins
>> 362
548 G4LinInterpolation linearAlgo; 363 G4LinInterpolation linearAlgo;
549 if (bin == 0) value = linearAlgo.Calculate(x 364 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
550 else value = algorithm->Calculate(x, bin, *p 365 else value = algorithm->Calculate(x, bin, *pdf, *energies);
551 366
>> 367 // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
552 return value; 368 return value;
553 } 369 }
554 370
555 //....oooOO0OOooo........oooOO0OOooo........oo <<
556 <<
557 G4double G4EMDataSet::IntegrationFunction(G4do 371 G4double G4EMDataSet::IntegrationFunction(G4double x)
558 { 372 {
559 // This function is needed by G4Integrator t 373 // This function is needed by G4Integrator to calculate the integral of the data distribution
>> 374
560 G4double y = 0; 375 G4double y = 0;
561 376
562 // Locate the random value in the X vector b << 377 // Locate the random value in the X vector based on the PDF
563 G4int bin = (G4int)FindLowerBound(x); << 378 size_t bin = FindLowerBound(x);
564 379
565 // Interpolate to calculate the X value: 380 // Interpolate to calculate the X value:
566 // linear interpolation in the first bin (to 381 // linear interpolation in the first bin (to avoid problem with 0),
567 // interpolation with associated algorithm i 382 // interpolation with associated algorithm in other bins
568 383
569 G4LinInterpolation linearAlgo; 384 G4LinInterpolation linearAlgo;
570 385
571 if (bin == 0) y = linearAlgo.Calculate(x, bi 386 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
572 else y = algorithm->Calculate(x, bin, *energ 387 else y = algorithm->Calculate(x, bin, *energies, *data);
573 388
574 return y; 389 return y;
575 } 390 }
576 391
577 392