Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4CrossSectionDataSet.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 64 
 65 G4CrossSectionDataSet::G4CrossSectionDataSet(G4VDataSetAlgorithm* argAlgorithm, 
 66                G4double argUnitEnergies, 
 67                G4double argUnitData) :
 68    algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
 69 {;}
 70 
 71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 72 
 73 G4CrossSectionDataSet::~G4CrossSectionDataSet()
 74 {
 75   CleanUpComponents();
 76  
 77   if (algorithm)
 78     delete algorithm;
 79 }
 80 
 81 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 82 
 83 G4bool G4CrossSectionDataSet::LoadData(const G4String & argFileName)
 84 {
 85   CleanUpComponents();
 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   G4int maxI = (G4int)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   G4int i(1);
203   while (i<maxI)
204     {
205       std::size_t 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       std::size_t 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 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
251 G4bool G4CrossSectionDataSet::LoadNonLogData(const G4String & argFileName)
252 {
253   CleanUpComponents();
254 
255   G4String fullFileName(FullFileName(argFileName));
256   std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
257 
258   if (!in.is_open())
259     {
260       G4String message("data file \"");
261       message+=fullFileName;
262       message+="\" not found";
263       G4Exception("G4CrossSectionDataSet::LoadNonLogData",
264         "em0003",FatalException,message);
265       return false;
266     }
267 
268   std::vector<G4DataVector *> columns;
269 
270   std::stringstream *stream(new std::stringstream);
271   char c;
272   G4bool comment(false);
273   G4bool space(true);
274   G4bool first(true);
275 
276   try
277     {
278       while (!in.eof())
279   {
280     in.get(c);
281    
282     switch (c)
283       {
284       case '\r':
285       case '\n':
286         if (!first)
287     {
288       unsigned long i(0);
289       G4double value;
290       
291       while (!stream->eof())
292         {
293           (*stream) >> value;
294        
295           while (i>=columns.size())
296                         {
297        columns.push_back(new G4DataVector);
298                         }
299       
300           columns[i]->push_back(value);
301        
302           i++;
303         }
304       
305       delete stream;
306       stream=new std::stringstream;
307     }
308      
309         first=true;
310         comment=false;
311         space=true;
312         break;
313      
314       case '#':
315         comment=true;
316         break;
317      
318       case '\t':
319       case ' ':
320         space = true;
321               break;
322 
323       default:
324         if (comment) { break; }
325         if (space && (!first)) { (*stream) << ' '; }
326 
327         first=false;
328         (*stream) << c;
329         space=false;
330       }
331   }
332     }
333   catch(const std::ios::failure &e)
334     {
335       // some implementations of STL could throw a "failture" exception
336       // when read wants read characters after end of file
337     }
338  
339   delete stream;
340  
341   G4int maxI = (G4int)columns.size();
342  
343   if (maxI<2)
344     {
345       G4String message("data file \"");
346       message+=fullFileName;
347       message+="\" should have at least two columns";
348       G4Exception("G4CrossSectionDataSet::LoadNonLogData",
349         "em0005",FatalException,message);
350       return false;
351     }
352  
353   G4int i(1);
354   while (i<maxI)
355     {
356       std::size_t maxJ(columns[i]->size());
357 
358       if (maxJ!=columns[0]->size())
359   {
360     G4String message("data file \"");
361     message+=fullFileName;
362     message+="\" has lines with a different number of columns";
363           G4Exception("G4CrossSectionDataSet::LoadNonLogData",
364         "em0005",FatalException,message);
365     return false;
366   }
367 
368       std::size_t j(0);
369 
370       G4DataVector *argEnergies=new G4DataVector;
371       G4DataVector *argData=new G4DataVector;
372 
373       while(j<maxJ)
374   {
375     argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
376     argData->push_back(columns[i]->operator[] (j)*GetUnitData());
377     j++;
378   }
379 
380       AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
381   
382       i++;
383     }
384 
385   i=maxI;
386   while (i>0)
387     {
388       i--;
389       delete columns[i];
390     }
391 
392   return true;
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
396 
397 G4bool G4CrossSectionDataSet::SaveData(const G4String & argFileName) const
398 {
399   const G4int n = (G4int)NumberOfComponents();
400  
401   if (n==0)
402     {
403       G4Exception("G4CrossSectionDataSet::SaveData",
404         "em0005",FatalException,"expected at least one component");
405       return false;
406     }
407  
408   G4String fullFileName(FullFileName(argFileName));
409   std::ofstream out(fullFileName);
410 
411   if (!out.is_open())
412     {
413       G4String message("cannot open \"");
414       message+=fullFileName;
415       message+="\"";
416       G4Exception("G4CrossSectionDataSet::SaveData",
417         "em0003",FatalException,message);
418       return false;
419     }
420  
421   G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
422   G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
423   G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
424  
425   G4int k(n);
426  
427   while (k>0)
428     {
429       k--;
430       iData[k]=GetComponent(k)->GetData(0).cbegin();
431     }
432  
433   while (iEnergies!=iEnergiesEnd)
434     {
435       out.precision(10);
436       out.width(15);
437       out.setf(std::ofstream::left);
438       out << ((*iEnergies)/GetUnitEnergies());
439   
440       k=0;
441   
442       while (k<n)
443   {
444     out << ' ';
445     out.precision(10);
446     out.width(15);
447     out.setf(std::ofstream::left);
448     out << ((*(iData[k]))/GetUnitData());
449 
450     iData[k]++;
451     k++;
452   }
453   
454       out << std::endl;
455   
456       iEnergies++;
457     }
458  
459   delete[] iData;
460 
461   return true;
462 }
463 
464 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
465 
466 G4String G4CrossSectionDataSet::FullFileName(const G4String& argFileName) const
467 {
468   const char* path = G4FindDataDir("G4LEDATA");
469   if (!path)
470     {
471       G4Exception("G4CrossSectionDataSet::FullFileName",
472        "em0006",FatalException,"G4LEDATA environment variable not set");
473       return "NULL";
474     }  
475 
476   std::ostringstream fullFileName;
477  
478   fullFileName << path << "/" << argFileName << ".dat";
479                       
480   return G4String(fullFileName.str().c_str());
481 }
482 
483 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
484 
485 G4double G4CrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
486 {
487   // Returns the sum over the shells corresponding to e
488   G4double value = 0.;
489 
490   std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
491   std::vector<G4VEMDataSet *>::const_iterator end(components.end());
492 
493   while (i!=end)
494     {
495       value+=(*i)->FindValue(argEnergy);
496       i++;
497     }
498 
499   return value;
500 }
501 
502 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
503 
504 void G4CrossSectionDataSet::PrintData(void) const
505 {
506   const G4int n = (G4int)NumberOfComponents();
507 
508   G4cout << "The data set has " << n << " components" << G4endl;
509   G4cout << G4endl;
510  
511   G4int i(0);
512  
513   while (i<n)
514     {
515       G4cout << "--- Component " << i << " ---" << G4endl;
516       GetComponent(i)->PrintData();
517       ++i;
518     }
519 }
520 
521 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
522 
523 void G4CrossSectionDataSet::SetEnergiesData(G4DataVector* argEnergies, 
524            G4DataVector* argData, 
525            G4int argComponentId)
526 {
527   G4VEMDataSet * component(components[argComponentId]);
528  
529   if (component)
530     {
531       component->SetEnergiesData(argEnergies, argData, 0);
532       return;
533     }
534 
535   std::ostringstream message;
536   message << "component " << argComponentId << " not found";
537  
538   G4Exception("G4CrossSectionDataSet::SetEnergiesData",
539        "em0005",FatalException,message.str().c_str());
540 }
541 
542 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
543 
544 void G4CrossSectionDataSet::SetLogEnergiesData(G4DataVector* argEnergies, 
545            G4DataVector* argData,
546                                          G4DataVector* argLogEnergies,
547                                          G4DataVector* argLogData, 
548            G4int argComponentId)
549 {
550   G4VEMDataSet * component(components[argComponentId]);
551  
552   if (component)
553     {
554       component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
555       return;
556     }
557 
558   std::ostringstream message;
559   message << "component " << argComponentId << " not found";
560  
561   G4Exception("G4CrossSectionDataSet::SetLogEnergiesData",
562        "em0005",FatalException,message.str().c_str());
563 }
564 
565 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
566 
567 void G4CrossSectionDataSet::CleanUpComponents()
568 {
569   while (!components.empty())
570     {
571       if (components.back()) delete components.back();
572       components.pop_back();
573     }
574 }
575 
576 
577