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