ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MicroElecCrossSectionDataSet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MicroElecCrossSectionDataSet.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 // 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 {
53 
54  if (algorithm)
55  delete algorithm;
56 }
57 
59 {
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("G4MicroElecCrossSectionDataSet::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  case ' ':
144  space = true;
145  break;
146 
147  default:
148  if (comment) { break; }
149  if (space && (!first)) { (*stream) << ' '; }
150 
151  first=false;
152  (*stream) << c;
153  space=false;
154  }
155  }
156  }
157  catch(const std::ios::failure &e)
158  {
159  // some implementations of STL could throw a "failture" exception
160  // when read wants read characters after end of file
161  }
162 
163  delete stream;
164 
165  std::vector<G4DataVector *>::size_type maxI(columns.size());
166 
167  if (maxI<2)
168  {
169  G4String message("Data file \"");
170  message+=fullFileName;
171  message+="\" should have at least two columns";
172  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
173  FatalException,message);
174  return false;
175  }
176 
177  std::vector<G4DataVector*>::size_type i(1);
178  while (i<maxI)
179  {
180  G4DataVector::size_type maxJ(columns[i]->size());
181 
182  if (maxJ!=columns[0]->size())
183  {
184  G4String message("Data file \"");
185  message+=fullFileName;
186  message+="\" has lines with a different number of columns";
187  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
188  FatalException,message);
189  return false;
190  }
191 
192  G4DataVector::size_type j(0);
193 
194  G4DataVector *argEnergies=new G4DataVector;
195  G4DataVector *argData=new G4DataVector;
196  G4DataVector *argLogEnergies=new G4DataVector;
197  G4DataVector *argLogData=new G4DataVector;
198 
199  while(j<maxJ)
200  {
201  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
202  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
203  argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
204  argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
205  j++;
206  }
207 
208  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
209 
210  i++;
211  }
212 
213  i=maxI;
214  while (i>0)
215  {
216  i--;
217  delete columns[i];
218  delete log_columns[i];
219  }
220 
221  return true;
222 }
223 
224 
226 {
228 
229  G4String fullFileName(FullFileName(argFileName));
230  std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
231 
232  if (!in.is_open())
233  {
234  G4String message("Data file \"");
235  message+=fullFileName;
236  message+="\" not found";
237  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0003",
238  FatalException,message);
239  return false;
240  }
241 
242  std::vector<G4DataVector *> columns;
243 
244  std::stringstream *stream(new std::stringstream);
245  char c;
246  G4bool comment(false);
247  G4bool space(true);
248  G4bool first(true);
249 
250  try
251  {
252  while (!in.eof())
253  {
254  in.get(c);
255 
256  switch (c)
257  {
258  case '\r':
259  case '\n':
260  if (!first)
261  {
262  unsigned long i(0);
263  G4double value;
264 
265  while (!stream->eof())
266  {
267  (*stream) >> value;
268 
269  while (i>=columns.size())
270  {
271  columns.push_back(new G4DataVector);
272  }
273 
274  columns[i]->push_back(value);
275 
276  i++;
277  }
278 
279  delete stream;
280  stream=new std::stringstream;
281  }
282 
283  first=true;
284  comment=false;
285  space=true;
286  break;
287 
288  case '#':
289  comment=true;
290  break;
291 
292  case '\t':
293  case ' ':
294  space = true;
295  break;
296 
297  default:
298  if (comment) { break; }
299  if (space && (!first)) { (*stream) << ' '; }
300 
301  first=false;
302  (*stream) << c;
303  space=false;
304  }
305  }
306  }
307  catch(const std::ios::failure &e)
308  {
309  // some implementations of STL could throw a "failture" exception
310  // when read wants read characters after end of file
311  }
312 
313  delete stream;
314 
315  std::vector<G4DataVector *>::size_type maxI(columns.size());
316 
317  if (maxI<2)
318  {
319  G4String message("Data file \"");
320  message+=fullFileName;
321  message+="\" should have at least two columns";
322  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
323  FatalException,message);
324  return false;
325  }
326 
327  std::vector<G4DataVector*>::size_type i(1);
328  while (i<maxI)
329  {
330  G4DataVector::size_type maxJ(columns[i]->size());
331 
332  if (maxJ!=columns[0]->size())
333  {
334  G4String message("Data file \"");
335  message+=fullFileName;
336  message+="\" has lines with a different number of columns.";
337  G4Exception("G4MicroElecCrossSectionDataSet::LoadData","em0005",
338  FatalException,message);
339  return false;
340  }
341 
342  G4DataVector::size_type j(0);
343 
344  G4DataVector *argEnergies=new G4DataVector;
345  G4DataVector *argData=new G4DataVector;
346 
347  while(j<maxJ)
348  {
349  argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
350  argData->push_back(columns[i]->operator[] (j)*GetUnitData());
351  j++;
352  }
353 
354  AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
355 
356  i++;
357  }
358 
359  i=maxI;
360  while (i>0)
361  {
362  i--;
363  delete columns[i];
364  }
365 
366  return true;
367 }
368 
369 
371 {
372  const size_t n(NumberOfComponents());
373 
374  if (n==0)
375  {
376  G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
377  FatalException,"Expected at least one component");
378 
379  return false;
380  }
381 
382  G4String fullFileName(FullFileName(argFileName));
383  std::ofstream out(fullFileName);
384 
385  if (!out.is_open())
386  {
387  G4String message("Cannot open \"");
388  message+=fullFileName;
389  message+="\"";
390  G4Exception("G4MicroElecCrossSectionDataSet::SaveData","em0005",
391  FatalException,message);
392  return false;
393  }
394 
395  G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
396  G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
397  G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
398 
399  size_t k(n);
400 
401  while (k>0)
402  {
403  k--;
404  iData[k]=GetComponent(k)->GetData(0).begin();
405  }
406 
407  while (iEnergies!=iEnergiesEnd)
408  {
409  out.precision(10);
410  out.width(15);
411  out.setf(std::ofstream::left);
412  out << ((*iEnergies)/GetUnitEnergies());
413 
414  k=0;
415 
416  while (k<n)
417  {
418  out << ' ';
419  out.precision(10);
420  out.width(15);
421  out.setf(std::ofstream::left);
422  out << ((*(iData[k]))/GetUnitData());
423 
424  iData[k]++;
425  k++;
426  }
427 
428  out << std::endl;
429 
430  iEnergies++;
431  }
432 
433  delete[] iData;
434 
435  return true;
436 }
437 
438 
440 {
441  char* path = std::getenv("G4LEDATA");
442  if (!path)
443  {
444  G4Exception("G4MicroElecCrossSectionDataSet::FullFileName","em0006",
445  FatalException,"G4LEDATA environment variable not set.");
446 
447  return "";
448  }
449 
450  std::ostringstream fullFileName;
451 
452  fullFileName << path << "/" << argFileName << ".dat";
453 
454  return G4String(fullFileName.str().c_str());
455 }
456 
457 
458 G4double G4MicroElecCrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
459 {
460  // Returns the sum over the shells corresponding to e
461  G4double value = 0.;
462 
463  std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
464  std::vector<G4VEMDataSet *>::const_iterator end(components.end());
465 
466  while (i!=end)
467  {
468  value+=(*i)->FindValue(argEnergy);
469  i++;
470  }
471 
472  return value;
473 }
474 
475 
477 {
478  const size_t n(NumberOfComponents());
479 
480  G4cout << "The data set has " << n << " components" << G4endl;
481  G4cout << G4endl;
482 
483  size_t i(0);
484 
485  while (i<n)
486  {
487  G4cout << "--- Component " << i << " ---" << G4endl;
488  GetComponent(i)->PrintData();
489  i++;
490  }
491 }
492 
493 
495  G4DataVector* argData,
496  G4int argComponentId)
497 {
498  G4VEMDataSet * component(components[argComponentId]);
499 
500  if (component)
501  {
502  component->SetEnergiesData(argEnergies, argData, 0);
503  return;
504  }
505 
506  std::ostringstream message;
507  message << "Component " << argComponentId << " not found";
508 
509  G4Exception("G4MicroElecCrossSectionDataSet::SetEnergiesData","em0005",
510  FatalException,message.str().c_str());
511 
512 }
513 
514 
516  G4DataVector* argData,
517  G4DataVector* argLogEnergies,
518  G4DataVector* argLogData,
519  G4int argComponentId)
520 {
521  G4VEMDataSet * component(components[argComponentId]);
522 
523  if (component)
524  {
525  component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
526  return;
527  }
528 
529  std::ostringstream message;
530  message << "Component " << argComponentId << " not found";
531 
532  G4Exception("G4MicroElecCrossSectionDataSet::SetLogEnergiesData","em0005",
533  FatalException,message.str().c_str());
534 
535 }
536 
537 
539 {
540  while (!components.empty())
541  {
542  if (components.back()) delete components.back();
543  components.pop_back();
544  }
545 }
546 
547