G4MuElecCrossSectionDataSet.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //  MR - 04/04/2012
00027 //  Based on G4DNACrossSectionDataSet
00028 // 
00029 
00030 
00031 #include "G4MuElecCrossSectionDataSet.hh"
00032 #include "G4VDataSetAlgorithm.hh"
00033 #include "G4EMDataSet.hh"
00034 #include <vector>
00035 #include <fstream>
00036 #include <sstream>
00037 
00038 
00039 G4MuElecCrossSectionDataSet::G4MuElecCrossSectionDataSet(G4VDataSetAlgorithm* argAlgorithm, 
00040                                                    G4double argUnitEnergies, 
00041                                                    G4double argUnitData)
00042   :
00043    algorithm(argAlgorithm), unitEnergies(argUnitEnergies), unitData(argUnitData)
00044 {
00045   z = 0;
00046 
00047 }
00048 
00049 
00050 G4MuElecCrossSectionDataSet::~G4MuElecCrossSectionDataSet()
00051 {
00052   CleanUpComponents();
00053  
00054   if (algorithm)
00055     delete algorithm;
00056 }
00057 
00058 G4bool G4MuElecCrossSectionDataSet::LoadData(const G4String & argFileName)
00059 {
00060   CleanUpComponents();
00061 
00062   G4String fullFileName(FullFileName(argFileName));
00063   std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
00064 
00065   if (!in.is_open())
00066     {
00067       G4String message("Data file \"");
00068       message+=fullFileName;
00069       message+="\" not found";
00070       G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0003",
00071                       FatalException,message);
00072       return false;
00073     }
00074 
00075   std::vector<G4DataVector *> columns;
00076   std::vector<G4DataVector *> log_columns;
00077 
00078   std::stringstream *stream(new std::stringstream);
00079   char c;
00080   G4bool comment(false);
00081   G4bool space(true);
00082   G4bool first(true);
00083 
00084   try
00085     {
00086       while (!in.eof())
00087         {
00088           in.get(c);
00089    
00090           switch (c)
00091             {
00092             case '\r':
00093             case '\n':
00094               if (!first)
00095                 {
00096                   unsigned long i(0);
00097                   G4double value;
00098       
00099                   while (!stream->eof())
00100                     {
00101                       (*stream) >> value;
00102        
00103                       while (i>=columns.size())
00104                         {
00105                          columns.push_back(new G4DataVector);
00106                          log_columns.push_back(new G4DataVector);
00107                         }
00108       
00109                       columns[i]->push_back(value);
00110 
00111 // N. A. Karakatsanis
00112 // A condition is applied to check if negative or zero values are present in the dataset.
00113 // If yes, then a near-zero value is applied to allow the computation of the logarithmic value
00114 // If a value is zero, this simplification is acceptable
00115 // If a value is negative, then it is not acceptable and the data of the particular column of
00116 // logarithmic values should not be used by interpolation methods.
00117 //
00118 // Therefore, G4LogLogInterpolation and G4LinLogLogInterpolation should not be used if negative values are present.
00119 // Instead, G4LinInterpolation is safe in every case
00120 // SemiLogInterpolation is safe only if the energy columns are non-negative
00121 // G4LinLogInterpolation is safe only if the cross section data columns are non-negative
00122 
00123                       if (value <=0.) value = 1e-300;
00124                       log_columns[i]->push_back(std::log10(value));
00125        
00126                       i++;
00127                     }
00128       
00129                   delete stream;
00130                   stream=new std::stringstream;
00131                 }
00132      
00133               first=true;
00134               comment=false;
00135               space=true;
00136               break;
00137      
00138             case '#':
00139               comment=true;
00140               break;
00141      
00142             case '\t':
00143               c=' ';
00144             case ' ':
00145               if (space)
00146                 break;
00147             default:
00148               if (comment)
00149                 break;
00150      
00151               if (c==' ')
00152                 space=true;
00153               else
00154                 {
00155                   if (space && (!first))
00156                     (*stream) << ' ';
00157       
00158                   first=false;
00159                   (*stream) << c;
00160                   space=false;
00161                 }
00162             }
00163         }
00164     }
00165   catch(const std::ios::failure &e)
00166     {
00167       // some implementations of STL could throw a "failture" exception
00168       // when read wants read characters after end of file
00169     }
00170  
00171   delete stream;
00172  
00173   std::vector<G4DataVector *>::size_type maxI(columns.size());
00174  
00175   if (maxI<2)
00176     {
00177       G4String message("Data file \"");
00178       message+=fullFileName;
00179       message+="\" should have at least two columns";
00180       G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0005",
00181                       FatalException,message);
00182       return false;
00183     }
00184  
00185   std::vector<G4DataVector*>::size_type i(1);
00186   while (i<maxI)
00187     {
00188       G4DataVector::size_type maxJ(columns[i]->size());
00189 
00190       if (maxJ!=columns[0]->size())
00191         {
00192           G4String message("Data file \"");
00193           message+=fullFileName;
00194           message+="\" has lines with a different number of columns";
00195           G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0005",
00196                       FatalException,message);
00197           return false;
00198         }
00199 
00200       G4DataVector::size_type j(0);
00201 
00202       G4DataVector *argEnergies=new G4DataVector;
00203       G4DataVector *argData=new G4DataVector;
00204       G4DataVector *argLogEnergies=new G4DataVector;
00205       G4DataVector *argLogData=new G4DataVector;
00206 
00207       while(j<maxJ)
00208         {
00209           argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
00210           argData->push_back(columns[i]->operator[] (j)*GetUnitData());
00211           argLogEnergies->push_back(log_columns[0]->operator[] (j) + std::log10(GetUnitEnergies()));
00212           argLogData->push_back(log_columns[i]->operator[] (j) + std::log10(GetUnitData()));
00213           j++;
00214         }
00215 
00216       AddComponent(new G4EMDataSet(i-1, argEnergies, argData, argLogEnergies, argLogData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
00217   
00218       i++;
00219     }
00220 
00221   i=maxI;
00222   while (i>0)
00223     {
00224       i--;
00225       delete columns[i];
00226       delete log_columns[i];
00227     }
00228 
00229   return true;
00230 }
00231 
00232 
00233 G4bool G4MuElecCrossSectionDataSet::LoadNonLogData(const G4String & argFileName)
00234 {
00235   CleanUpComponents();
00236 
00237   G4String fullFileName(FullFileName(argFileName));
00238   std::ifstream in(fullFileName, std::ifstream::binary|std::ifstream::in);
00239 
00240   if (!in.is_open())
00241     {
00242       G4String message("Data file \"");
00243       message+=fullFileName;
00244       message+="\" not found";
00245       G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0003",
00246                       FatalException,message);
00247       return false;
00248     }
00249 
00250   std::vector<G4DataVector *> columns;
00251 
00252   std::stringstream *stream(new std::stringstream);
00253   char c;
00254   G4bool comment(false);
00255   G4bool space(true);
00256   G4bool first(true);
00257 
00258   try
00259     {
00260       while (!in.eof())
00261         {
00262           in.get(c);
00263    
00264           switch (c)
00265             {
00266             case '\r':
00267             case '\n':
00268               if (!first)
00269                 {
00270                   unsigned long i(0);
00271                   G4double value;
00272       
00273                   while (!stream->eof())
00274                     {
00275                       (*stream) >> value;
00276        
00277                       while (i>=columns.size())
00278                         {
00279                          columns.push_back(new G4DataVector);
00280                         }
00281       
00282                       columns[i]->push_back(value);
00283        
00284                       i++;
00285                     }
00286       
00287                   delete stream;
00288                   stream=new std::stringstream;
00289                 }
00290      
00291               first=true;
00292               comment=false;
00293               space=true;
00294               break;
00295      
00296             case '#':
00297               comment=true;
00298               break;
00299      
00300             case '\t':
00301               c=' ';
00302             case ' ':
00303               if (space)
00304                 break;
00305             default:
00306               if (comment)
00307                 break;
00308      
00309               if (c==' ')
00310                 space=true;
00311               else
00312                 {
00313                   if (space && (!first))
00314                     (*stream) << ' ';
00315       
00316                   first=false;
00317                   (*stream) << c;
00318                   space=false;
00319                 }
00320             }
00321         }
00322     }
00323   catch(const std::ios::failure &e)
00324     {
00325       // some implementations of STL could throw a "failture" exception
00326       // when read wants read characters after end of file
00327     }
00328  
00329   delete stream;
00330  
00331   std::vector<G4DataVector *>::size_type maxI(columns.size());
00332  
00333   if (maxI<2)
00334     {
00335       G4String message("Data file \"");
00336       message+=fullFileName;
00337       message+="\" should have at least two columns";
00338       G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0005",
00339                       FatalException,message);
00340       return false;
00341     }
00342  
00343   std::vector<G4DataVector*>::size_type i(1);
00344   while (i<maxI)
00345     {
00346       G4DataVector::size_type maxJ(columns[i]->size());
00347 
00348       if (maxJ!=columns[0]->size())
00349         {
00350           G4String message("Data file \"");
00351           message+=fullFileName;
00352           message+="\" has lines with a different number of columns.";
00353           G4Exception("G4MuElecCrossSectionDataSet::LoadData","em0005",
00354                       FatalException,message);
00355           return false;
00356         }
00357 
00358       G4DataVector::size_type j(0);
00359 
00360       G4DataVector *argEnergies=new G4DataVector;
00361       G4DataVector *argData=new G4DataVector;
00362 
00363       while(j<maxJ)
00364         {
00365           argEnergies->push_back(columns[0]->operator[] (j)*GetUnitEnergies());
00366           argData->push_back(columns[i]->operator[] (j)*GetUnitData());
00367           j++;
00368         }
00369 
00370       AddComponent(new G4EMDataSet(i-1, argEnergies, argData, GetAlgorithm()->Clone(), GetUnitEnergies(), GetUnitData()));
00371   
00372       i++;
00373     }
00374 
00375   i=maxI;
00376   while (i>0)
00377     {
00378       i--;
00379       delete columns[i];
00380     }
00381 
00382   return true;
00383 }
00384 
00385 
00386 G4bool G4MuElecCrossSectionDataSet::SaveData(const G4String & argFileName) const
00387 {
00388   const size_t n(NumberOfComponents());
00389  
00390   if (n==0)
00391     {
00392       G4Exception("G4MuElecCrossSectionDataSet::SaveData","em0005",
00393                       FatalException,"Expected at least one component");
00394 
00395       return false;
00396     }
00397  
00398   G4String fullFileName(FullFileName(argFileName));
00399   std::ofstream out(fullFileName);
00400 
00401   if (!out.is_open())
00402     {
00403       G4String message("Cannot open \"");
00404       message+=fullFileName;
00405       message+="\"";
00406       G4Exception("G4MuElecCrossSectionDataSet::SaveData","em0005",
00407                       FatalException,message);
00408       return false;
00409     }
00410  
00411   G4DataVector::const_iterator iEnergies(GetComponent(0)->GetEnergies(0).begin());
00412   G4DataVector::const_iterator iEnergiesEnd(GetComponent(0)->GetEnergies(0).end());
00413   G4DataVector::const_iterator * iData(new G4DataVector::const_iterator[n]);
00414  
00415   size_t k(n);
00416  
00417   while (k>0)
00418     {
00419       k--;
00420       iData[k]=GetComponent(k)->GetData(0).begin();
00421     }
00422  
00423   while (iEnergies!=iEnergiesEnd)
00424     {
00425       out.precision(10);
00426       out.width(15);
00427       out.setf(std::ofstream::left);
00428       out << ((*iEnergies)/GetUnitEnergies());
00429   
00430       k=0;
00431   
00432       while (k<n)
00433         {
00434           out << ' ';
00435           out.precision(10);
00436           out.width(15);
00437           out.setf(std::ofstream::left);
00438           out << ((*(iData[k]))/GetUnitData());
00439 
00440           iData[k]++;
00441           k++;
00442         }
00443   
00444       out << std::endl;
00445   
00446       iEnergies++;
00447     }
00448  
00449   delete[] iData;
00450 
00451   return true;
00452 }
00453 
00454 
00455 G4String G4MuElecCrossSectionDataSet::FullFileName(const G4String& argFileName) const
00456 {
00457   char* path = getenv("G4LEDATA");
00458   if (!path)
00459   {
00460       G4Exception("G4MuElecCrossSectionDataSet::FullFileName","em0006",
00461                       FatalException,"G4LEDATA environment variable not set.");
00462       
00463       return "";
00464   }
00465   
00466   std::ostringstream fullFileName;
00467  
00468   fullFileName << path << "/" << argFileName << ".dat";
00469                       
00470   return G4String(fullFileName.str().c_str());
00471 }
00472 
00473 
00474 G4double G4MuElecCrossSectionDataSet::FindValue(G4double argEnergy, G4int /* argComponentId */) const
00475 {
00476   // Returns the sum over the shells corresponding to e
00477   G4double value = 0.;
00478 
00479   std::vector<G4VEMDataSet *>::const_iterator i(components.begin());
00480   std::vector<G4VEMDataSet *>::const_iterator end(components.end());
00481 
00482   while (i!=end)
00483     {
00484       value+=(*i)->FindValue(argEnergy);
00485       i++;
00486     }
00487 
00488   return value;
00489 }
00490 
00491 
00492 void G4MuElecCrossSectionDataSet::PrintData(void) const
00493 {
00494   const size_t n(NumberOfComponents());
00495 
00496   G4cout << "The data set has " << n << " components" << G4endl;
00497   G4cout << G4endl;
00498  
00499   size_t i(0);
00500  
00501   while (i<n)
00502     {
00503       G4cout << "--- Component " << i << " ---" << G4endl;
00504       GetComponent(i)->PrintData();
00505       i++;
00506     }
00507 }
00508 
00509 
00510 void G4MuElecCrossSectionDataSet::SetEnergiesData(G4DataVector* argEnergies, 
00511                                          G4DataVector* argData, 
00512                                          G4int argComponentId)
00513 {
00514   G4VEMDataSet * component(components[argComponentId]);
00515  
00516   if (component)
00517     {
00518       component->SetEnergiesData(argEnergies, argData, 0);
00519       return;
00520     }
00521 
00522   std::ostringstream message;
00523   message << "Component " << argComponentId << " not found";
00524  
00525   G4Exception("G4MuElecCrossSectionDataSet::SetEnergiesData","em0005",
00526                  FatalException,message.str().c_str());
00527   
00528 }
00529 
00530 
00531 void G4MuElecCrossSectionDataSet::SetLogEnergiesData(G4DataVector* argEnergies, 
00532                                          G4DataVector* argData,
00533                                          G4DataVector* argLogEnergies,
00534                                          G4DataVector* argLogData, 
00535                                          G4int argComponentId)
00536 {
00537   G4VEMDataSet * component(components[argComponentId]);
00538  
00539   if (component)
00540     {
00541       component->SetLogEnergiesData(argEnergies, argData, argLogEnergies, argLogData, 0);
00542       return;
00543     }
00544 
00545   std::ostringstream message;
00546   message << "Component " << argComponentId << " not found";
00547  
00548   G4Exception("G4MuElecCrossSectionDataSet::SetLogEnergiesData","em0005",
00549                  FatalException,message.str().c_str());
00550 
00551 }
00552 
00553 
00554 void G4MuElecCrossSectionDataSet::CleanUpComponents()
00555 {
00556   while (!components.empty())
00557     {
00558       if (components.back()) delete components.back();
00559       components.pop_back();
00560     }
00561 }
00562 
00563 

Generated on Mon May 27 17:48:53 2013 for Geant4 by  doxygen 1.4.7