G4EMDataSet.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 // History:
00027 // -----------
00028 // 31 Jul 2001   MGP        Created
00029 //
00030 // 15 Jul 2009   Nicolas A. Karakatsanis
00031 //
00032 //                           - New Constructor was added when logarithmic data are loaded as well
00033 //                             to enhance CPU performance
00034 //
00035 //                           - Destructor was modified accordingly
00036 //
00037 //                           - LoadNonLogData method was created to load only the non-logarithmic data from G4EMLOW
00038 //                             dataset. It is essentially performing the data loading operations as in the past.
00039 //
00040 //                           - LoadData method was revised in order to calculate the logarithmic values of the data
00041 //                             It retrieves the data values from the G4EMLOW data files but, then, calculates the
00042 //                             respective log values and loads them to seperate data structures. 
00043 //
00044 //                           - SetLogEnergiesData method was cretaed to set logarithmic values to G4 data vectors.
00045 //                             The EM data sets, initialized this way, contain both non-log and log values.
00046 //                             These initialized data sets can enhance the computing performance of data interpolation
00047 //                             operations
00048 //
00049 // 26 Dec 2010 V.Ivanchenko Fixed Coverity warnings and cleanup logic
00050 //
00051 // -------------------------------------------------------------------
00052 
00053 #include "G4EMDataSet.hh"
00054 #include "G4VDataSetAlgorithm.hh"
00055 #include <fstream>
00056 #include <sstream>
00057 #include "G4Integrator.hh"
00058 #include "Randomize.hh"
00059 #include "G4LinInterpolation.hh"
00060 
00061 
00062 G4EMDataSet::G4EMDataSet(G4int Z, 
00063                          G4VDataSetAlgorithm* algo, 
00064                          G4double xUnit, 
00065                          G4double yUnit,
00066                          G4bool random): 
00067   z(Z),
00068   energies(0),
00069   data(0),
00070   log_energies(0),
00071   log_data(0),
00072   algorithm(algo),
00073   unitEnergies(xUnit),
00074   unitData(yUnit),
00075   pdf(0),
00076   randomSet(random)
00077 {
00078   if (algorithm == 0) {
00079     G4Exception("G4EMDataSet::G4EMDataSet",
00080                     "em1012",FatalException,"interpolation == 0");
00081   } else if (randomSet) { BuildPdf(); }
00082 }
00083 
00084 G4EMDataSet::G4EMDataSet(G4int argZ, 
00085                          G4DataVector* dataX, 
00086                          G4DataVector* dataY, 
00087                          G4VDataSetAlgorithm* algo, 
00088                          G4double xUnit, 
00089                          G4double yUnit,
00090                          G4bool random):
00091   z(argZ),
00092   energies(dataX),
00093   data(dataY),
00094   log_energies(0),
00095   log_data(0),
00096   algorithm(algo),
00097   unitEnergies(xUnit),
00098   unitData(yUnit),
00099   pdf(0),
00100   randomSet(random)
00101 {
00102   if (algorithm == 0) {
00103     G4Exception("G4EMDataSet::G4EMDataSet",
00104                     "em1012",FatalException,"interpolation == 0");
00105   } else if ((energies == 0) ^ (data == 0)) {
00106     G4Exception("G4EMDataSet::G4EMDataSet",
00107                     "em1012",FatalException,"different size for energies and data (zero case)");
00108   } else if (energies->size() != data->size()) { 
00109     G4Exception("G4EMDataSet::G4EMDataSet",
00110                     "em1012",FatalException,"different size for energies and data");
00111   } else if (randomSet) {
00112     BuildPdf();
00113   }
00114 }
00115 
00116 G4EMDataSet::G4EMDataSet(G4int argZ, 
00117                          G4DataVector* dataX, 
00118                          G4DataVector* dataY,
00119                          G4DataVector* dataLogX,
00120                          G4DataVector* dataLogY, 
00121                          G4VDataSetAlgorithm* algo, 
00122                          G4double xUnit, 
00123                          G4double yUnit,
00124                          G4bool random):
00125   z(argZ),
00126   energies(dataX),
00127   data(dataY),
00128   log_energies(dataLogX),
00129   log_data(dataLogY),
00130   algorithm(algo),
00131   unitEnergies(xUnit),
00132   unitData(yUnit),
00133   pdf(0),
00134   randomSet(random)
00135 {
00136   if (algorithm == 0) {
00137     G4Exception("G4EMDataSet::G4EMDataSet",
00138                     "em1012",FatalException,"interpolation == 0");
00139   } else if ((energies == 0) ^ (data == 0)) {
00140     G4Exception("G4EMDataSet::G4EMDataSet",
00141                     "em1012",FatalException,"different size for energies and data (zero case)");
00142   } else if (energies->size() != data->size()) { 
00143     G4Exception("G4EMDataSet::G4EMDataSet",
00144                     "em1012",FatalException,"different size for energies and data");
00145   } else if ((log_energies == 0) ^ (log_data == 0)) {
00146     G4Exception("G4EMDataSet::G4EMDataSet",
00147                     "em1012",FatalException,"different size for log energies and log data (zero case)");
00148   } else if (log_energies->size() != log_data->size()) {
00149     G4Exception("G4EMDataSet::G4EMDataSet",
00150                     "em1012",FatalException,"different size for log energies and log data");
00151   } else if (randomSet) {
00152     BuildPdf();
00153   }
00154 }
00155 
00156 
00157 G4EMDataSet::~G4EMDataSet()
00158 { 
00159   delete algorithm;
00160   delete energies; 
00161   delete data; 
00162   delete pdf; 
00163   delete log_energies; 
00164   delete log_data; 
00165 }
00166 
00167 
00168 G4double G4EMDataSet::FindValue(G4double energy, G4int /* componentId */) const
00169 {
00170   if (!energies) {
00171    G4Exception("G4EMDataSet::FindValue",
00172                     "em1012",FatalException,"energies == 0");
00173     return 0.0;
00174   } else if (energies->empty()) {
00175     return 0.0;
00176   } else if (energy <= (*energies)[0]) {
00177     return (*data)[0];
00178   }
00179   size_t i = energies->size()-1;
00180   if (energy >= (*energies)[i]) { return (*data)[i]; }
00181 
00182   //Nicolas A. Karakatsanis: Check if the logarithmic data have been loaded to decide
00183   //                         which Interpolation-Calculation method will be applied
00184   if (log_energies != 0) 
00185    {
00186      return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data, *log_energies, *log_data);
00187    }
00188 
00189   return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
00190 }
00191 
00192 
00193 void G4EMDataSet::PrintData(void) const
00194 {
00195   if (!energies)
00196     {
00197       G4cout << "Data not available." << G4endl;
00198     }
00199   else
00200     {
00201       size_t size = energies->size();
00202       for (size_t i(0); i<size; i++)
00203         {
00204           G4cout << "Point: " << ((*energies)[i]/unitEnergies)
00205                  << " - Data value: " << ((*data)[i]/unitData);
00206           if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
00207           G4cout << G4endl; 
00208         }
00209     }
00210 }
00211 
00212 
00213 void G4EMDataSet::SetEnergiesData(G4DataVector* dataX, 
00214                                   G4DataVector* dataY, 
00215                                   G4int /* componentId */)
00216 {
00217   if (energies) { delete energies; }
00218   energies = dataX;
00219 
00220   if (data) { delete data; }
00221   data = dataY;
00222  
00223   if ((energies == 0) ^ (data==0)) {
00224     G4Exception("G4EMDataSet::SetEnergiesData",
00225                     "em1012",FatalException,"different size for energies and data (zero case)");
00226     return;
00227   } else if (energies == 0) { return; }
00228 
00229   //G4cout << "Size of energies: " << energies->size() << G4endl << "Size of data: " << data->size() << G4endl;
00230   if (energies->size() != data->size()) { 
00231    G4Exception("G4EMDataSet::SetEnergiesData",
00232                     "em1012",FatalException,"different size for energies and data");
00233   }
00234 }
00235 
00236 void G4EMDataSet::SetLogEnergiesData(G4DataVector* dataX,
00237                                      G4DataVector* dataY,
00238                                      G4DataVector* data_logX, 
00239                                      G4DataVector* data_logY,
00240                                      G4int /* componentId */)
00241 {
00242   //Load of the actual energy and data values  
00243   if (energies) { delete energies; }
00244   energies = dataX;
00245   if (data) { delete data; }
00246   data = dataY;
00247   //Load of the logarithmic energy and data values
00248   if (log_energies) { delete log_energies; }
00249   log_energies = data_logX;
00250   if (log_data) { delete log_data; }
00251   log_data = data_logY;
00252 
00253   //Check if data loaded properly from data files
00254   if ( !energies ) {
00255     if(data || log_energies || log_data ) { 
00256       G4Exception("G4EMDataSet::SetLogEnergiesData",
00257                     "em1012",FatalException,"inconsistent data");
00258     }
00259     return;
00260   } else {
00261     if ( !data ) {
00262       G4Exception("G4EMDataSet::SetLogEnergiesData",
00263                     "em1012",FatalException,"only energy, no data");
00264       return; 
00265     } else if (energies->size() != data->size()) { 
00266       G4Exception("G4EMDataSet::SetLogEnergiesData",
00267                     "em1012",FatalException,"different size for energies and data");
00268       return;
00269     }
00270     //G4cout << "Size of energies: " << energies->size() << G4endl << "Size of data: " << data->size() << G4endl << G4endl;
00271 
00272     //Check if logarithmic data loaded properly from data files 
00273     if ( !log_energies ) {
00274       if(log_data) {
00275         G4Exception("G4EMDataSet::SetLogEnergiesData",
00276                     "em1012",FatalException,"inconsistence of log_data");
00277       }
00278       return;
00279     } else {
00280       if ( !log_data ) { 
00281         G4Exception("G4EMDataSet::SetLogEnergiesData",
00282                     "em1012",FatalException,"only log_energies, no data");
00283       } else if ((log_energies->size() != log_data->size()) || (log_energies->size() != data->size())) { 
00284         G4Exception("G4EMDataSet::SetLogEnergiesData",
00285                     "em1012",FatalException,"different size for log energies and data");
00286       }
00287     }
00288   }
00289   //G4cout << "Size of log energies: " << log_energies->size() << G4endl << "Size of log data: " << log_data->size() << G4endl;  
00290 }
00291 
00292 G4bool G4EMDataSet::LoadData(const G4String& fileName)
00293 {
00294  // The file is organized into four columns:
00295  // 1st column contains the values of energy
00296  // 2nd column contains the corresponding data value
00297  // The file terminates with the pattern: -1   -1
00298  //                                       -2   -2
00299 
00300   G4String fullFileName(FullFileName(fileName));
00301   std::ifstream in(fullFileName);
00302 
00303   if (!in.is_open())
00304     {
00305       G4String message("data file \"");
00306       message += fullFileName;
00307       message += "\" not found";
00308       G4Exception("G4EMDataSet::LoadData",
00309                     "em1012",FatalException,message);
00310       return false;
00311     }
00312 
00313   delete energies;
00314   delete data;   
00315   delete log_energies;
00316   delete log_data;   
00317   energies = new G4DataVector;
00318   data = new G4DataVector;
00319   log_energies = new G4DataVector;
00320   log_data = new G4DataVector;
00321 
00322   G4double a, b;
00323   //G4int k = 0;
00324   //G4int nColumns = 2;
00325 
00326   do
00327     {
00328       in >> a >> b;
00329   
00330       if (a != -1 && a != -2)
00331         {
00332           if (a==0.) { a=1e-300; }
00333           if (b==0.) { b=1e-300; }
00334           a *= unitEnergies;
00335           b *= unitData;
00336           energies->push_back(a);
00337           log_energies->push_back(std::log10(a));
00338           data->push_back(b);
00339           log_data->push_back(std::log10(b));
00340         }
00341     }
00342   while (a != -2);
00343 
00344   if (randomSet) { BuildPdf(); }
00345  
00346   return true;
00347 }
00348 
00349 
00350 G4bool G4EMDataSet::LoadNonLogData(const G4String& fileName)
00351 {
00352  // The file is organized into four columns:
00353  // 1st column contains the values of energy
00354  // 2nd column contains the corresponding data value
00355  // The file terminates with the pattern: -1   -1
00356  //                                       -2   -2
00357 
00358   G4String fullFileName(FullFileName(fileName));
00359   std::ifstream in(fullFileName);
00360   if (!in.is_open())
00361     {
00362       G4String message("data file \"");
00363       message += fullFileName;
00364       message += "\" not found";
00365       G4Exception("G4EMDataSet::LoadNonLogData",
00366                     "em1012",FatalException,message);
00367     }
00368 
00369   G4DataVector* argEnergies=new G4DataVector;
00370   G4DataVector* argData=new G4DataVector;
00371 
00372   G4double a;
00373   G4int k = 0;
00374   G4int nColumns = 2;
00375 
00376   do
00377     {
00378       in >> a;
00379   
00380       if (a != -1 && a != -2)
00381         {
00382           if (k%nColumns == 0)
00383             {
00384              argEnergies->push_back(a*unitEnergies);
00385             }
00386           else if (k%nColumns == 1)
00387             {
00388              argData->push_back(a*unitData);
00389             }
00390           k++;
00391         }
00392     }
00393   while (a != -2);
00394  
00395   SetEnergiesData(argEnergies, argData, 0);
00396 
00397   if (randomSet) BuildPdf();
00398  
00399   return true;
00400 }
00401 
00402 
00403 
00404 G4bool G4EMDataSet::SaveData(const G4String& name) const
00405 {
00406   // The file is organized into two columns:
00407   // 1st column is the energy
00408   // 2nd column is the corresponding value
00409   // The file terminates with the pattern: -1   -1
00410   //                                       -2   -2
00411  
00412   G4String fullFileName(FullFileName(name));
00413   std::ofstream out(fullFileName);
00414 
00415   if (!out.is_open())
00416     {
00417       G4String message("cannot open \"");
00418       message+=fullFileName;
00419       message+="\"";
00420       G4Exception("G4EMDataSet::SaveData",
00421                     "em1012",FatalException,message);
00422     }
00423  
00424   out.precision(10);
00425   out.width(15);
00426   out.setf(std::ofstream::left);
00427   
00428   if (energies!=0 && data!=0)
00429     {
00430       G4DataVector::const_iterator i(energies->begin());
00431       G4DataVector::const_iterator endI(energies->end());
00432       G4DataVector::const_iterator j(data->begin());
00433   
00434       while (i!=endI)
00435         {
00436           out.precision(10);
00437           out.width(15);
00438           out.setf(std::ofstream::left);
00439           out << ((*i)/unitEnergies) << ' ';
00440 
00441           out.precision(10);
00442           out.width(15);
00443           out.setf(std::ofstream::left);
00444           out << ((*j)/unitData) << std::endl;
00445 
00446           i++;
00447           j++;
00448         }
00449     }
00450  
00451   out.precision(10);
00452   out.width(15);
00453   out.setf(std::ofstream::left);
00454   out << -1.f << ' ';
00455 
00456   out.precision(10);
00457   out.width(15);
00458   out.setf(std::ofstream::left);
00459   out << -1.f << std::endl;
00460 
00461   out.precision(10);
00462   out.width(15);
00463   out.setf(std::ofstream::left);
00464   out << -2.f << ' ';
00465 
00466   out.precision(10);
00467   out.width(15);
00468   out.setf(std::ofstream::left);
00469   out << -2.f << std::endl;
00470 
00471   return true;
00472 }
00473 
00474 
00475 
00476 size_t G4EMDataSet::FindLowerBound(G4double x) const
00477 {
00478   size_t lowerBound = 0;
00479   size_t upperBound(energies->size() - 1);
00480   
00481   while (lowerBound <= upperBound) 
00482     {
00483       size_t midBin((lowerBound + upperBound) / 2);
00484 
00485       if (x < (*energies)[midBin]) upperBound = midBin - 1;
00486       else lowerBound = midBin + 1;
00487     }
00488   
00489   return upperBound;
00490 }
00491 
00492 
00493 size_t G4EMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
00494 {
00495   size_t lowerBound = 0;;
00496   size_t upperBound(values->size() - 1);
00497   
00498   while (lowerBound <= upperBound) 
00499     {
00500       size_t midBin((lowerBound + upperBound) / 2);
00501 
00502       if (x < (*values)[midBin]) upperBound = midBin - 1;
00503       else lowerBound = midBin + 1;
00504     }
00505   
00506   return upperBound;
00507 }
00508 
00509 
00510 G4String G4EMDataSet::FullFileName(const G4String& name) const
00511 {
00512   char* path = getenv("G4LEDATA");
00513   if (!path) {
00514      G4Exception("G4EMDataSet::FullFileName",
00515                     "em0006",FatalException,"G4LEDATA environment variable not set");
00516     return "";
00517   }
00518   std::ostringstream fullFileName;
00519   fullFileName << path << '/' << name << z << ".dat";
00520                       
00521   return G4String(fullFileName.str().c_str());
00522 }
00523 
00524 
00525 
00526 void G4EMDataSet::BuildPdf() 
00527 {
00528   pdf = new G4DataVector;
00529   G4Integrator <G4EMDataSet, G4double(G4EMDataSet::*)(G4double)> integrator;
00530 
00531   G4int nData = data->size();
00532   pdf->push_back(0.);
00533 
00534   // Integrate the data distribution 
00535   G4int i;
00536   G4double totalSum = 0.;
00537   for (i=1; i<nData; i++)
00538     {
00539       G4double xLow = (*energies)[i-1];
00540       G4double xHigh = (*energies)[i];
00541       G4double sum = integrator.Legendre96(this, &G4EMDataSet::IntegrationFunction, xLow, xHigh);
00542       totalSum = totalSum + sum;
00543       pdf->push_back(totalSum);
00544     }
00545 
00546   // Normalize to the last bin
00547   G4double tot = 0.;
00548   if (totalSum > 0.) tot = 1. / totalSum;
00549   for (i=1;  i<nData; i++)
00550     {
00551       (*pdf)[i] = (*pdf)[i] * tot;
00552     }
00553 }
00554 
00555 G4double G4EMDataSet::RandomSelect(G4int /* componentId */) const
00556 {
00557   G4double value = 0.;
00558   // Random select a X value according to the cumulative probability distribution
00559   // derived from the data
00560 
00561   if (!pdf) {
00562     G4Exception("G4EMDataSet::RandomSelect",
00563                     "em1012",FatalException,"PDF has not been created for this data set");
00564     return value;
00565   }
00566 
00567   G4double x = G4UniformRand();
00568 
00569   // Locate the random value in the X vector based on the PDF
00570   size_t bin = FindLowerBound(x,pdf);
00571 
00572   // Interpolate the PDF to calculate the X value: 
00573   // linear interpolation in the first bin (to avoid problem with 0),
00574   // interpolation with associated data set algorithm in other bins
00575 
00576   G4LinInterpolation linearAlgo;
00577   if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
00578   else value = algorithm->Calculate(x, bin, *pdf, *energies);
00579 
00580   //  G4cout << x << " random bin "<< bin << " - " << value << G4endl;
00581   return value;
00582 }
00583 
00584 G4double G4EMDataSet::IntegrationFunction(G4double x)
00585 {
00586   // This function is needed by G4Integrator to calculate the integral of the data distribution
00587 
00588   G4double y = 0;
00589 
00590  // Locate the random value in the X vector based on the PDF
00591   size_t bin = FindLowerBound(x);
00592 
00593   // Interpolate to calculate the X value: 
00594   // linear interpolation in the first bin (to avoid problem with 0),
00595   // interpolation with associated algorithm in other bins
00596 
00597   G4LinInterpolation linearAlgo;
00598   
00599   if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
00600   else y = algorithm->Calculate(x, bin, *energies, *data);
00601  
00602   return y;
00603 }
00604 

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