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