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