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