Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 105122 2017-07-13 13:47:51Z gcosmo $
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 {
78  CleanUpComponents();
79 
80  if (algorithm)
81  delete algorithm;
82 }
83 
85 {
86  CleanUpComponents();
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  case ' ':
170  space = true;
171  break;
172 
173  default:
174  if (comment) { break; }
175  if (space && (!first)) { (*stream) << ' '; }
176 
177  first=false;
178  (*stream) << c;
179  space=false;
180  }
181  }
182  }
183  catch(const std::ios::failure &e)
184  {
185  // some implementations of STL could throw a "failture" exception
186  // when read wants read characters after end of file
187  }
188 
189  delete stream;
190 
191  std::vector<G4DataVector *>::size_type maxI(columns.size());
192 
193  if (maxI<2)
194  {
195  G4String message("data file \"");
196  message+=fullFileName;
197  message+="\" should have at least two columns";
198  G4Exception("G4CrossSectionDataSet::LoadData",
199  "em0005",FatalException,message);
200  return false;
201  }
202 
203  std::vector<G4DataVector*>::size_type i(1);
204  while (i<maxI)
205  {
206  G4DataVector::size_type maxJ(columns[i]->size());
207 
208  if (maxJ!=columns[0]->size())
209  {
210  G4String message("data file \"");
211  message+=fullFileName;
212  message+="\" has lines with a different number of columns";
213  G4Exception("G4CrossSectionDataSet::LoadData",
214  "em0005",FatalException,message);
215  return false;
216  }
217 
218  G4DataVector::size_type j(0);
219 
220  G4DataVector *argEnergies=new G4DataVector;
221  G4DataVector *argData=new G4DataVector;
222  G4DataVector *argLogEnergies=new G4DataVector;
223  G4DataVector *argLogData=new G4DataVector;
224 
225  while(j<maxJ)
226  {
227  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
228  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
229  argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
230  argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
231  j++;
232  }
233 
234  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
235 
236  i++;
237  }
238 
239  i=maxI;
240  while (i>0)
241  {
242  i--;
243  delete columns[i];
244  delete log_columns[i];
245  }
246 
247  return true;
248 }
249 
250 
252 {
253  CleanUpComponents();
254 
255  G4String fullFileName(FullFileName(argFileName));
256  std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
257 
258  if (!in.is_open())
259  {
260  G4String message("data file \"");
261  message+=fullFileName;
262  message+="\" not found";
263  G4Exception("G4CrossSectionDataSet::LoadNonLogData",
264  "em0003",FatalException,message);
265  return false;
266  }
267 
268  std::vector<G4DataVector *> columns;
269 
270  std::stringstream *stream(new std::stringstream);
271  char c;
272  G4bool comment(false);
273  G4bool space(true);
274  G4bool first(true);
275 
276  try
277  {
278  while (!in.eof())
279  {
280  in.get(c);
281 
282  switch (c)
283  {
284  case '\r':
285  case '\n':
286  if (!first)
287  {
288  unsigned long i(0);
289  G4double value;
290 
291  while (!stream->eof())
292  {
293  (*stream) >> value;
294 
295  while (i>=columns.size())
296  {
297  columns.push_back(new G4DataVector);
298  }
299 
300  columns[i]->push_back(value);
301 
302  i++;
303  }
304 
305  delete stream;
306  stream=new std::stringstream;
307  }
308 
309  first=true;
310  comment=false;
311  space=true;
312  break;
313 
314  case '#':
315  comment=true;
316  break;
317 
318  case '\t':
319  case ' ':
320  space = true;
321  break;
322 
323  default:
324  if (comment) { break; }
325  if (space && (!first)) { (*stream) << ' '; }
326 
327  first=false;
328  (*stream) << c;
329  space=false;
330  }
331  }
332  }
333  catch(const std::ios::failure &e)
334  {
335  // some implementations of STL could throw a "failture" exception
336  // when read wants read characters after end of file
337  }
338 
339  delete stream;
340 
341  std::vector<G4DataVector *>::size_type maxI(columns.size());
342 
343  if (maxI<2)
344  {
345  G4String message("data file \"");
346  message+=fullFileName;
347  message+="\" should have at least two columns";
348  G4Exception("G4CrossSectionDataSet::LoadNonLogData",
349  "em0005",FatalException,message);
350  return false;
351  }
352 
353  std::vector<G4DataVector*>::size_type i(1);
354  while (i<maxI)
355  {
356  G4DataVector::size_type maxJ(columns[i]->size());
357 
358  if (maxJ!=columns[0]->size())
359  {
360  G4String message("data file \"");
361  message+=fullFileName;
362  message+="\" has lines with a different number of columns";
363  G4Exception("G4CrossSectionDataSet::LoadNonLogData",
364  "em0005",FatalException,message);
365  return false;
366  }
367 
368  G4DataVector::size_type j(0);
369 
370  G4DataVector *argEnergies=new G4DataVector;
371  G4DataVector *argData=new G4DataVector;
372 
373  while(j<maxJ)
374  {
375  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
376  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
377  j++;
378  }
379 
380  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
381 
382  i++;
383  }
384 
385  i=maxI;
386  while (i>0)
387  {
388  i--;
389  delete columns[i];
390  }
391 
392  return true;
393 }
394 
395 
397 {
398  const size_t n(NumberOfComponents());
399 
400  if (n==0)
401  {
402  G4Exception("G4CrossSectionDataSet::SaveData",
403  "em0005",FatalException,"expected at least one component");
404  return false;
405  }
406 
407  G4String fullFileName(FullFileName(argFileName));
408  std::ofstream out(fullFileName);
409 
410  if (!out.is_open())
411  {
412  G4String message("cannot open \"");
413  message+=fullFileName;
414  message+="\"";
415  G4Exception("G4CrossSectionDataSet::SaveData",
416  "em0003",FatalException,message);
417  return false;
418  }
419 
420  G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
421  G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
422  G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
423 
424  size_t k(n);
425 
426  while (k>0)
427  {
428  k--;
429  iData[k]=GetComponent(k)->GetData(0).begin();
430  }
431 
432  while (iEnergies!=iEnergiesEnd)
433  {
434  out.precision(10);
435  out.width(15);
436  out.setf(std::ofstream::left);
437  out << ((*iEnergies)/GetUnitEnergies());
438 
439  k=0;
440 
441  while (k<n)
442  {
443  out << ' ';
444  out.precision(10);
445  out.width(15);
446  out.setf(std::ofstream::left);
447  out << ((*(iData[k]))/GetUnitData());
448 
449  iData[k]++;
450  k++;
451  }
452 
453  out << std::endl;
454 
455  iEnergies++;
456  }
457 
458  delete[] iData;
459 
460  return true;
461 }
462 
463 
464 G4String G4CrossSectionDataSet::FullFileName(const G4String& argFileName) const
465 {
466  char* path = getenv("G4LEDATA");
467  if (!path)
468  {
469  G4Exception("G4CrossSectionDataSet::FullFileName",
470  "em0006",FatalException,"G4LEDATA environment variable not set");
471  return "NULL";
472  }
473 
474  std::ostringstream fullFileName;
475 
476  fullFileName << path << "/" << argFileName << ".dat";
477 
478  return G4String(fullFileName.str().c_str());
479 }
480 
481 
482 G4double G4CrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
483 {
484  // Returns the sum over the shells corresponding to e
485  G4double value = 0.;
486 
487  std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
488  std::vector<G4VEMDataSet *>::const_iterator end(components.end());
489 
490  while (i!=end)
491  {
492  value+=(*i)->FindValue(argEnergy);
493  i++;
494  }
495 
496  return value;
497 }
498 
499 
501 {
502  const size_t n(NumberOfComponents());
503 
504  G4cout << "The data set has " << n << " components" << G4endl;
505  G4cout << G4endl;
506 
507  size_t i(0);
508 
509  while (i<n)
510  {
511  G4cout << "--- Component " << i << " ---" << G4endl;
512  GetComponent(i)->PrintData();
513  i++;
514  }
515 }
516 
517 
519  G4DataVector* argData,
520  G4int argComponentId)
521 {
522  G4VEMDataSet * component(components[argComponentId]);
523 
524  if (component)
525  {
526  component->SetEnergiesData(argEnergies, argData, 0);
527  return;
528  }
529 
530  std::ostringstream message;
531  message << "component " << argComponentId << " not found";
532 
533  G4Exception("G4CrossSectionDataSet::SetEnergiesData",
534  "em0005",FatalException,message.str().c_str());
535 }
536 
537 
539  G4DataVector* argData,
540  G4DataVector* argLogEnergies,
541  G4DataVector* argLogData,
542  G4int argComponentId)
543 {
544  G4VEMDataSet * component(components[argComponentId]);
545 
546  if (component)
547  {
548  component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
549  return;
550  }
551 
552  std::ostringstream message;
553  message << "component " << argComponentId << " not found";
554 
555  G4Exception("G4CrossSectionDataSet::SetLogEnergiesData",
556  "em0005",FatalException,message.str().c_str());
557 }
558 
559 
560 void G4CrossSectionDataSet::CleanUpComponents()
561 {
562  while (!components.empty())
563  {
564  if (components.back()) delete components.back();
565  components.pop_back();
566  }
567 }
568 
569 
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
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
virtual void PrintData(void) const
G4CrossSectionDataSet(G4VDataSetAlgorithm *algo, G4double xUnit=CLHEP::MeV, G4double dataUnit=CLHEP::barn)
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
bool G4bool
Definition: G4Types.hh:79
virtual G4double FindValue(G4double e, G4int componentId=0) 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
virtual const G4VEMDataSet * GetComponent(G4int componentId) const
virtual void AddComponent(G4VEMDataSet *dataSet)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
virtual void SetLogEnergiesData(G4DataVector *x, G4DataVector *values, G4DataVector *log_x, G4DataVector *log_values, G4int componentId)
virtual G4bool LoadData(const G4String &argFileName)