Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MicroElecCrossSectionDataSet.cc
Go to the documentation of this file.
1 
2 //
3 // ********************************************************************
4 // * License and Disclaimer *
5 // * *
6 // * The Geant4 software is copyright of the Copyright Holders of *
7 // * the Geant4 Collaboration. It is provided under the terms and *
8 // * conditions of the Geant4 Software License, included in the file *
9 // * LICENSE and available at http://cern.ch/geant4/license . These *
10 // * include a list of copyright holders. *
11 // * *
12 // * Neither the authors of this software system, nor their employing *
13 // * institutes,nor the agencies providing financial support for this *
14 // * work make any representation or warranty, express or implied, *
15 // * regarding this software system or assume any liability for its *
16 // * use. Please see the license in the file LICENSE and URL above *
17 // * for the full disclaimer and the limitation of liability. *
18 // * *
19 // * This code implementation is the result of the scientific and *
20 // * technical work of the GEANT4 collaboration. *
21 // * By using, copying, modifying or distributing the software (or *
22 // * any work based on the software) you agree to acknowledge its *
23 // * use in resulting scientific publications, and indicate your *
24 // * acceptance of all terms of the Geant4 Software license. *
25 // ********************************************************************
26 //
27 // MR - 04/04/2012
28 // Based on G4DNACrossSectionDataSet
29 //
30 
31 
33 #include "G4VDataSetAlgorithm.hh"
34 #include "G4EMDataSet.hh"
35 #include <vector>
36 #include <fstream>
37 #include <sstream>
38 
39 
41  G4double argUnitEnergies,
42  G4double argUnitData)
43  :
44  algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
45 {
46  z = 0;
47 
48 }
49 
50 
52 {
53  CleanUpComponents();
54 
55  if (algorithm)
56  delete algorithm;
57 }
58 
60 {
61  CleanUpComponents();
62 
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::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);
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  case ' ':
145  space = true;
146  break;
147 
148  default:
149  if (comment) { break; }
150  if (space && (!first)) { (*stream) << ' '; }
151 
152  first=false;
153  (*stream) << c;
154  space=false;
155  }
156  }
157  }
158  catch(const std::ios::failure &e)
159  {
160  // some implementations of STL could throw a "failture" exception
161  // when read wants read characters after end of file
162  }
163 
164  delete stream;
165 
166  std::vector<G4DataVector *>::size_type maxI(columns.size());
167 
168  if (maxI<2)
169  {
170  G4String message("Data file \"");
171  message+=fullFileName;
172  message+="\" should have at least two columns";
173  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
174  FatalException,message);
175  return false;
176  }
177 
178  std::vector<G4DataVector*>::size_type i(1);
179  while (i<maxI)
180  {
181  G4DataVector::size_type maxJ(columns[i]->size());
182 
183  if (maxJ!=columns[0]->size())
184  {
185  G4String message("Data file \"");
186  message+=fullFileName;
187  message+="\" has lines with a different number of columns";
188  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
189  FatalException,message);
190  return false;
191  }
192 
193  G4DataVector::size_type j(0);
194 
195  G4DataVector *argEnergies=new G4DataVector;
196  G4DataVector *argData=new G4DataVector;
197  G4DataVector *argLogEnergies=new G4DataVector;
198  G4DataVector *argLogData=new G4DataVector;
199 
200  while(j<maxJ)
201  {
202  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
203  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
204  argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
205  argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
206  j++;
207  }
208 
209  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
210 
211  i++;
212  }
213 
214  i=maxI;
215  while (i>0)
216  {
217  i--;
218  delete columns[i];
219  delete log_columns[i];
220  }
221 
222  return true;
223 }
224 
225 
227 {
228  CleanUpComponents();
229 
230  G4String fullFileName(FullFileName(argFileName));
231  std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
232 
233  if (!in.is_open())
234  {
235  G4String message("Data file \"");
236  message+=fullFileName;
237  message+="\" not found";
238  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0003",
239  FatalException,message);
240  return false;
241  }
242 
243  std::vector<G4DataVector *> columns;
244 
245  std::stringstream *stream(new std::stringstream);
246  char c;
247  G4bool comment(false);
248  G4bool space(true);
249  G4bool first(true);
250 
251  try
252  {
253  while (!in.eof())
254  {
255  in.get(c);
256 
257  switch (c)
258  {
259  case '\r':
260  case '\n':
261  if (!first)
262  {
263  unsigned long i(0);
264  G4double value;
265 
266  while (!stream->eof())
267  {
268  (*stream) >> value;
269 
270  while (i>=columns.size())
271  {
272  columns.push_back(new G4DataVector);
273  }
274 
275  columns[i]->push_back(value);
276 
277  i++;
278  }
279 
280  delete stream;
281  stream=new std::stringstream;
282  }
283 
284  first=true;
285  comment=false;
286  space=true;
287  break;
288 
289  case '#':
290  comment=true;
291  break;
292 
293  case '\t':
294  case ' ':
295  space = true;
296  break;
297 
298  default:
299  if (comment) { break; }
300  if (space && (!first)) { (*stream) << ' '; }
301 
302  first=false;
303  (*stream) << c;
304  space=false;
305  }
306  }
307  }
308  catch(const std::ios::failure &e)
309  {
310  // some implementations of STL could throw a "failture" exception
311  // when read wants read characters after end of file
312  }
313 
314  delete stream;
315 
316  std::vector<G4DataVector *>::size_type maxI(columns.size());
317 
318  if (maxI<2)
319  {
320  G4String message("Data file \"");
321  message+=fullFileName;
322  message+="\" should have at least two columns";
323  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
324  FatalException,message);
325  return false;
326  }
327 
328  std::vector<G4DataVector*>::size_type i(1);
329  while (i<maxI)
330  {
331  G4DataVector::size_type maxJ(columns[i]->size());
332 
333  if (maxJ!=columns[0]->size())
334  {
335  G4String message("Data file \"");
336  message+=fullFileName;
337  message+="\" has lines with a different number of columns.";
338  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
339  FatalException,message);
340  return false;
341  }
342 
343  G4DataVector::size_type j(0);
344 
345  G4DataVector *argEnergies=new G4DataVector;
346  G4DataVector *argData=new G4DataVector;
347 
348  while(j<maxJ)
349  {
350  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
351  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
352  j++;
353  }
354 
355  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
356 
357  i++;
358  }
359 
360  i=maxI;
361  while (i>0)
362  {
363  i--;
364  delete columns[i];
365  }
366 
367  return true;
368 }
369 
370 
372 {
373  const size_t n(NumberOfComponents());
374 
375  if (n==0)
376  {
377  G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
378  FatalException,"Expected at least one component");
379 
380  return false;
381  }
382 
383  G4String fullFileName(FullFileName(argFileName));
384  std::ofstream out(fullFileName);
385 
386  if (!out.is_open())
387  {
388  G4String message("Cannot open \"");
389  message+=fullFileName;
390  message+="\"";
391  G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
392  FatalException,message);
393  return false;
394  }
395 
396  G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
397  G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
398  G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
399 
400  size_t k(n);
401 
402  while (k>0)
403  {
404  k--;
405  iData[k]=GetComponent(k)->GetData(0).begin();
406  }
407 
408  while (iEnergies!=iEnergiesEnd)
409  {
410  out.precision(10);
411  out.width(15);
412  out.setf(std::ofstream::left);
413  out << ((*iEnergies)/GetUnitEnergies());
414 
415  k=0;
416 
417  while (k<n)
418  {
419  out << ' ';
420  out.precision(10);
421  out.width(15);
422  out.setf(std::ofstream::left);
423  out << ((*(iData[k]))/GetUnitData());
424 
425  iData[k]++;
426  k++;
427  }
428 
429  out << std::endl;
430 
431  iEnergies++;
432  }
433 
434  delete[] iData;
435 
436  return true;
437 }
438 
439 
440 G4String G4MicroElecCrossSectionDataSet::FullFileName(const G4String& argFileName) const
441 {
442  char* path = getenv("G4LEDATA");
443  if (!path)
444  {
445  G4Exception("G4MicroElecCrossSectionDataSet::FullFileName","em0006",
446  FatalException,"G4LEDATA environment variable not set.");
447 
448  return "";
449  }
450 
451  std::ostringstream fullFileName;
452 
453  fullFileName << path << "/" << argFileName << ".dat";
454 
455  return G4String(fullFileName.str().c_str());
456 }
457 
458 
459 G4double G4MicroElecCrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
460 {
461  // Returns the sum over the shells corresponding to e
462  G4double value = 0.;
463 
464  std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
465  std::vector<G4VEMDataSet *>::const_iterator end(components.end());
466 
467  while (i!=end)
468  {
469  value+=(*i)->FindValue(argEnergy);
470  i++;
471  }
472 
473  return value;
474 }
475 
476 
478 {
479  const size_t n(NumberOfComponents());
480 
481  G4cout << "The data set has " << n << " components" << G4endl;
482  G4cout << G4endl;
483 
484  size_t i(0);
485 
486  while (i<n)
487  {
488  G4cout << "--- Component " << i << " ---" << G4endl;
489  GetComponent(i)->PrintData();
490  i++;
491  }
492 }
493 
494 
496  G4DataVector* argData,
497  G4int argComponentId)
498 {
499  G4VEMDataSet * component(components[argComponentId]);
500 
501  if (component)
502  {
503  component->SetEnergiesData(argEnergies, argData, 0);
504  return;
505  }
506 
507  std::ostringstream message;
508  message << "Component " << argComponentId << " not found";
509 
510  G4Exception("G4MicroElecCrossSectionDataSet::SetEnergiesData","em0005",
511  FatalException,message.str().c_str());
512 
513 }
514 
515 
517  G4DataVector* argData,
518  G4DataVector* argLogEnergies,
519  G4DataVector* argLogData,
520  G4int argComponentId)
521 {
522  G4VEMDataSet * component(components[argComponentId]);
523 
524  if (component)
525  {
526  component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
527  return;
528  }
529 
530  std::ostringstream message;
531  message << "Component " << argComponentId << " not found";
532 
533  G4Exception("G4MicroElecCrossSectionDataSet::SetLogEnergiesData","em0005",
534  FatalException,message.str().c_str());
535 
536 }
537 
538 
539 void G4MicroElecCrossSectionDataSet::CleanUpComponents()
540 {
541  while (!components.empty())
542  {
543  if (components.back()) delete components.back();
544  components.pop_back();
545  }
546 }
547 
548 
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual void SetEnergiesData(G4DataVector *x, G4DataVector *data, G4int component=0)=0
virtual void SetEnergiesData(G4DataVector *x, G4DataVector *values, G4int componentId)
virtual G4bool LoadNonLogData(const G4String &argFileName)
virtual G4bool SaveData(const G4String &argFileName) const
virtual const G4DataVector & GetData(G4int componentId) const =0
G4MicroElecCrossSectionDataSet(G4VDataSetAlgorithm *algo, G4double xUnit=CLHEP::MeV, G4double dataUnit=CLHEP::barn)
int G4int
Definition: G4Types.hh:78
virtual void SetLogEnergiesData(G4DataVector *x, G4DataVector *values, G4DataVector *log_x, G4DataVector *log_values, G4int componentId)
virtual void AddComponent(G4VEMDataSet *dataSet)
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
virtual size_t NumberOfComponents(void) const
bool G4bool
Definition: G4Types.hh:79
const G4int n
virtual void SetLogEnergiesData(G4DataVector *x, G4DataVector *data, G4DataVector *Log_x, G4DataVector *Log_data, G4int component=0)=0
virtual G4double FindValue(G4double e, G4int componentId=0) const
virtual G4bool LoadData(const G4String &argFileName)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
virtual void PrintData(void) const =0
virtual const G4DataVector & GetEnergies(G4int componentId) const
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13