ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuElecCrossSectionDataSet.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuElecCrossSectionDataSet.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 
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 {
49 
50  if (algorithm)
51  delete algorithm;
52 }
53 
55 {
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 {
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 
436 {
437  char* path = std::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 
535 {
536  while (!components.empty())
537  {
538  if (components.back()) delete components.back();
539  components.pop_back();
540  }
541 }
542 
543