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