| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // MR - 04/04/2012
27 // Based on G4DNACrossSectionDataSet
28 //
29
30
31 #include "G4MicroElecCrossSectionDataSet.hh"
32 #include "G4VDataSetAlgorithm.hh"
33 #include "G4EMDataSet.hh"
34 #include <vector>
35 #include <fstream>
36 #include <sstream>
37
38
39 G4MicroElecCrossSectionDataSet::G4MicroElecCrossSectionDataSet(G4VDataSetAlgorithm* argAlgorithm,
40 G4double argUnitEnergies,
41 G4double argUnitData)
42 :
43 algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
44 {;}
45
46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47 G4MicroElecCrossSectionDataSet::~G4MicroElecCrossSectionDataSet()
48 {
49 CleanUpComponents();
50
51 if (algorithm)
52 delete algorithm;
53 }
54
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 G4bool G4MicroElecCrossSectionDataSet::LoadData(const G4String & argFileName)
57 {
58 CleanUpComponents();
59
60 G4String fullFileName(FullFileName(argFileName));
61 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
62
63 if (!in.is_open())
64 {
65 G4String message("Data file \"");
66 message+=fullFileName;
67 message+="\" not found";
68 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0003",
69 FatalException,message);
70 return false;
71 }
72
73 std::vector<G4DataVector *> columns;
74 std::vector<G4DataVector *> log_columns;
75
76 std::stringstream *stream(new std::stringstream);
77 char c;
78 G4bool comment(false);
79 G4bool space(true);
80 G4bool first(true);
81
82 try
83 {
84 while (!in.eof())
85 {
86 in.get(c);
87
88 switch (c)
89 {
90 case '\r':
91 case '\n':
92 if (!first)
93 {
94 unsigned long i(0);
95 G4double value;
96
97 while (!stream->eof())
98 {
99 (*stream) >> value;
100
101 while (i>=columns.size())
102 {
103 columns.push_back(new G4DataVector);
104 log_columns.push_back(new G4DataVector);
105 }
106
107 columns[i]->push_back(value);
108 // N. A. Karakatsanis
109 // A condition is applied to check if negative or zero values are present in the dataset.
110 // If yes, then a near-zero value is applied to allow the computation of the logarithmic value
111 // If a value is zero, this simplification is acceptable
112 // If a value is negative, then it is not acceptable and the data of the particular column of
113 // logarithmic values should not be used by interpolation methods.
114 //
115 // Therefore, G4LogLogInterpolation and G4LinLogLogInterpolation should not be used if negative values are present.
116 // Instead, G4LinInterpolation is safe in every case
117 // SemiLogInterpolation is safe only if the energy columns are non-negative
118 // G4LinLogInterpolation is safe only if the cross section data columns are non-negative
119
120 if (value <=0.) value = 1e-300;
121 log_columns[i]->push_back(std::log10(value));
122
123 ++i;
124 }
125 delete stream;
126 stream=new std::stringstream;
127 }
128 first=true;
129 comment=false;
130 space=true;
131 break;
132
133 case '#':
134 comment=true;
135 break;
136 case '\t':
137 case ' ':
138 space = true;
139 break;
140 default:
141 if (comment) { break; }
142 if (space && (!first)) { (*stream) << ' '; }
143 first=false;
144 (*stream) << c;
145 space=false;
146 }
147 }
148 }
149 catch(const std::ios::failure &e)
150 {
151 // some implementations of STL could throw a "failture" exception
152 // when read wants read characters after end of file
153 }
154
155 delete stream;
156
157 std::size_t maxI(columns.size());
158
159 if (maxI<2)
160 {
161 G4String message("Data file \"");
162 message+=fullFileName;
163 message+="\" should have at least two columns";
164 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
165 FatalException,message);
166 return false;
167 }
168
169 std::size_t i(1);
170 while (i<maxI)
171 {
172 std::size_t maxJ(columns[i]->size());
173
174 if (maxJ!=columns[0]->size())
175 {
176 G4String message("Data file \"");
177 message+=fullFileName;
178 message+="\" has lines with a different number of columns";
179 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
180 FatalException,message);
181 return false;
182 }
183
184 std::size_t j(0);
185
186 G4DataVector *argEnergies=new G4DataVector;
187 G4DataVector *argData=new G4DataVector;
188 G4DataVector *argLogEnergies=new G4DataVector;
189 G4DataVector *argLogData=new G4DataVector;
190
191 while(j<maxJ)
192 {
193 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
194 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
195 argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
196 argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
197 j++;
198 }
199
200 AddComponent(new G4EMDataSet(G4int(i-1), argEnergies, argData, argLogEnergies, argLogData,
201 GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
202
203 ++i;
204 }
205
206 i=maxI;
207 while (i>0)
208 {
209 --i;
210 delete columns[i];
211 delete log_columns[i];
212 }
213
214 return true;
215 }
216
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218
219 G4bool G4MicroElecCrossSectionDataSet::LoadNonLogData(const G4String & argFileName)
220 {
221 CleanUpComponents();
222
223 G4String fullFileName(FullFileName(argFileName));
224 std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
225
226 if (!in.is_open())
227 {
228 G4String message("Data file \"");
229 message+=fullFileName;
230 message+="\" not found";
231 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0003",
232 FatalException,message);
233 return false;
234 }
235
236 std::vector<G4DataVector *> columns;
237
238 std::stringstream *stream(new std::stringstream);
239 char c;
240 G4bool comment(false);
241 G4bool space(true);
242 G4bool first(true);
243
244 try
245 {
246 while (!in.eof())
247 {
248 in.get(c);
249
250 switch (c)
251 {
252 case '\r':
253 case '\n':
254 if (!first)
255 {
256 unsigned long i(0);
257 G4double value;
258
259 while (!stream->eof())
260 {
261 (*stream) >> value;
262
263 while (i>=columns.size())
264 {
265 columns.push_back(new G4DataVector);
266 }
267
268 columns[i]->push_back(value);
269
270 i++;
271 }
272
273 delete stream;
274 stream=new std::stringstream;
275 }
276
277 first=true;
278 comment=false;
279 space=true;
280 break;
281
282 case '#':
283 comment=true;
284 break;
285
286 case '\t':
287 case ' ':
288 space = true;
289 break;
290
291 default:
292 if (comment) { break; }
293 if (space && (!first)) { (*stream) << ' '; }
294
295 first=false;
296 (*stream) << c;
297 space=false;
298 }
299 }
300 }
301 catch(const std::ios::failure &e)
302 {
303 // some implementations of STL could throw a "failture" exception
304 // when read wants read characters after end of file
305 }
306
307 delete stream;
308
309 std::size_t maxI(columns.size());
310
311 if (maxI<2)
312 {
313 G4String message("Data file \"");
314 message+=fullFileName;
315 message+="\" should have at least two columns";
316 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
317 FatalException,message);
318 return false;
319 }
320
321 std::size_t i(1);
322 while (i<maxI)
323 {
324 std::size_t maxJ(columns[i]->size());
325
326 if (maxJ!=columns[0]->size())
327 {
328 G4String message("Data file \"");
329 message+=fullFileName;
330 message+="\" has lines with a different number of columns.";
331 G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
332 FatalException,message);
333 return false;
334 }
335
336 std::size_t j(0);
337
338 G4DataVector *argEnergies=new G4DataVector;
339 G4DataVector *argData=new G4DataVector;
340
341 while(j<maxJ)
342 {
343 argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
344 argData->push_back(columns[i]->operator[] (j)*GetUnitData());
345 j++;
346 }
347
348 AddComponent(new G4EMDataSet(G4int(i-1), argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
349
350 i++;
351 }
352
353 i=maxI;
354 while (i>0)
355 {
356 i--;
357 delete columns[i];
358 }
359
360 return true;
361 }
362
363 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
364
365 G4bool G4MicroElecCrossSectionDataSet::SaveData(const G4String & argFileName) const
366 {
367 const std::size_t n(NumberOfComponents());
368 if (n==0)
369 {
370 G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
371 FatalException,"Expected at least one component");
372
373 return false;
374 }
375
376 G4String fullFileName(FullFileName(argFileName));
377 std::ofstream out(fullFileName);
378
379 if (!out.is_open())
380 {
381 G4String message("Cannot open \"");
382 message+=fullFileName;
383 message+="\"";
384 G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
385 FatalException,message);
386 return false;
387 }
388
389 G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
390 G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
391 G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
392
393 std::size_t k(n);
394
395 while (k>0)
396 {
397 k--;
398 iData[k]=GetComponent((G4int)k)->GetData(0).begin();
399 }
400
401 while (iEnergies!=iEnergiesEnd)
402 {
403 out.precision(10);
404 out.width(15);
405 out.setf(std::ofstream::left);
406 out << ((*iEnergies)/GetUnitEnergies());
407
408 k=0;
409
410 while (k<n)
411 {
412 out << ' ';
413 out.precision(10);
414 out.width(15);
415 out.setf(std::ofstream::left);
416 out << ((*(iData[k]))/GetUnitData());
417
418 iData[k]++;
419 k++;
420 }
421 out << std::endl;
422 iEnergies++;
423 }
424 delete[] iData;
425
426 return true;
427 }
428
429 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
430
431 G4String G4MicroElecCrossSectionDataSet::FullFileName(const G4String& argFileName) const
432 {
433 const char* path = G4FindDataDir("G4LEDATA");
434 if (!path)
435 {
436 G4Exception("G4MicroElecCrossSectionDataSet::FullFileName","em0006",
437 FatalException,"G4LEDATA environment variable not set.");
438
439 return "";
440 }
441
442 std::ostringstream fullFileName;
443
444 fullFileName << path << "/" << argFileName << ".dat";
445
446 return G4String(fullFileName.str().c_str());
447 }
448
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450
451 G4double G4MicroElecCrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
452 {
453 // Returns the sum over the shells corresponding to e
454 G4double value = 0.;
455
456 std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
457 std::vector<G4VEMDataSet *>::const_iterator end(components.end());
458
459 while (i!=end)
460 {
461 value+=(*i)->FindValue(argEnergy);
462 i++;
463 }
464
465 return value;
466 }
467
468 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
469
470 void G4MicroElecCrossSectionDataSet::PrintData(void) const
471 {
472 const std::size_t n(NumberOfComponents());
473
474 G4cout << "The data set has " << n << " components" << G4endl;
475 G4cout << G4endl;
476
477 G4int i(0);
478
479 while (i<(G4int)n)
480 {
481 G4cout << "--- Component " << i << " ---" << G4endl;
482 GetComponent(i)->PrintData();
483 i++;
484 }
485 }
486
487 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
488
489 void G4MicroElecCrossSectionDataSet::SetEnergiesData(G4DataVector* argEnergies,
490 G4DataVector* argData,
491 G4int argComponentId)
492 {
493 G4VEMDataSet * component(components[argComponentId]);
494
495 if (component)
496 {
497 component->SetEnergiesData(argEnergies, argData, 0);
498 return;
499 }
500
501 std::ostringstream message;
502 message << "Component " << argComponentId << " not found";
503
504 G4Exception("G4MicroElecCrossSectionDataSet::SetEnergiesData","em0005",
505 FatalException,message.str().c_str());
506
507 }
508
509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
510
511 void G4MicroElecCrossSectionDataSet::SetLogEnergiesData(G4DataVector* argEnergies,
512 G4DataVector* argData,
513 G4DataVector* argLogEnergies,
514 G4DataVector* argLogData,
515 G4int argComponentId)
516 {
517 G4VEMDataSet * component(components[argComponentId]);
518
519 if (component)
520 {
521 component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
522 return;
523 }
524
525 std::ostringstream message;
526 message << "Component " << argComponentId << " not found";
527
528 G4Exception("G4MicroElecCrossSectionDataSet::SetLogEnergiesData","em0005",
529 FatalException,message.str().c_str());
530
531 }
532
533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
534
535 void G4MicroElecCrossSectionDataSet::CleanUpComponents()
536 {
537 while (!components.empty())
538 {
539 if (components.back()) delete components.back();
540 components.pop_back();
541 }
542 }
543
544
545