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