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