G4LogLogInterpolation.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 //
00027 // $Id$
00028 //
00029 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
00030 //         Sebastian Incerti (incerti@cenbg.in2p3.fr)
00031 //         Nicolas A. Karakatsanis (knicolas@mail.ntua.gr)
00032 // History:
00033 // -----------
00034 // 31 Jul 2001   MGP        Created
00035 // 27 Jun 2008   SI         Add check to avoid FPE errors
00036 // 08 Dec 2008   NAK        Log-Log interpolation math formula streamlined, self-test function
00037 // 14 Jun 2008   NAK        New implementation for log-log interpolation after directly loading
00038 //                          logarithmic values from G4EMLOW dataset
00039 // -------------------------------------------------------------------
00040 
00041 #include "G4LogLogInterpolation.hh"
00042 
00043 // Constructor
00044 
00045 G4LogLogInterpolation::G4LogLogInterpolation()
00046 { }
00047 
00048 // Destructor
00049 
00050 G4LogLogInterpolation::~G4LogLogInterpolation()
00051 { }
00052 
00053 G4VDataSetAlgorithm* G4LogLogInterpolation::Clone() const 
00054 { return new G4LogLogInterpolation; }
00055 
00056 
00057 G4double G4LogLogInterpolation::Calculate(G4double x, G4int bin, 
00058                                           const G4DataVector& points, 
00059                                           const G4DataVector& data) const
00060 {
00061   //G4cout << "G4LogLogInterpolation is performed (2 arguments) " << G4endl;
00062   G4int nBins = data.size() - 1;
00063 //G4double oldresult = 0.;
00064   G4double value = 0.;
00065   if (x < points[0])
00066     {
00067       value = 0.;
00068     }
00069   else if (bin < nBins)
00070     {
00071       G4double e1 = points[bin];
00072       G4double e2 = points[bin+1];
00073       G4double d1 = data[bin];
00074       G4double d2 = data[bin+1];
00075 // Check of e1, e2, d1 and d2 values to avoid floating-point errors when estimating the interpolated value below -- S.I., Jun. 2008
00076       if ((d1 > 0.) && (d2 > 0.) && (e1 > 0.) && (e2 > 0.))
00077         {
00078 // Streamline the Log-Log Interpolation formula in order to reduce the required number of log10() function calls
00079 // Variable oldresult contains the result of old implementation of Log-Log interpolation -- M.G.P. Jun. 2001
00080 //       oldresult = (std::log10(d1)*std::log10(e2/x) + std::log10(d2)*std::log10(x/e1)) / std::log10(e2/e1);
00081 //       oldresult = std::pow(10.,oldresult);
00082 // Variable value contains the result of new implementation, after streamlining the math operation -- N.A.K. Oct. 2008
00083          value = std::log10(d1)+(std::log10(d2/d1)/std::log10(e2/e1)*std::log10(x/e1));
00084          value = std::pow(10.,value);
00085 // Test of the new implementation result (value variable) against the old one (oldresult) -- N.A.K. Dec. 2008
00086 //       G4double diffResult = value - oldresult;
00087 //       G4double relativeDiff = 1e-11;
00088 // Comparison of the two values based on a max allowable relative difference
00089 //       if ( std::fabs(diffResult) > relativeDiff*std::fabs(oldresult) )
00090 //        {
00091 // Abort comparison when at least one of two results is infinite
00092 //           if ((!std::isinf(oldresult)) && (!std::isinf(value)))
00093 //            {
00094 //              G4cout << "G4LogLogInterpolation> Old Interpolated Value is:" << oldresult << G4endl;
00095 //              G4cout << "G4LogLogInterpolation> New Interpolated Value is:" << value << G4endl << G4endl;
00096 //              G4cerr << "G4LogLogInterpolation> Error in Interpolation:" << G4endl;
00097 //              G4cerr << "The difference between new and old interpolated value is:" << diffResult << G4endl << G4endl;
00098 //            }
00099 //        }
00100         }
00101       else value = 0.;
00102     }
00103   else
00104     {
00105       value = data[nBins];
00106     }
00107   return value;
00108 }
00109 
00110 
00111 // Nicolas A. Karakatsanis: New implementation of log-log interpolation after directly loading 
00112 //                          logarithmic values from G4EMLOW dataset
00113 
00114 G4double G4LogLogInterpolation::Calculate(G4double x, G4int bin, 
00115                                           const G4DataVector& points,
00116                                           const G4DataVector& data,
00117                                           const G4DataVector& log_points, 
00118                                           const G4DataVector& log_data) const
00119 {
00120   //G4cout << "G4LogLogInterpolation is performed (4 arguments) " << G4endl;
00121   G4int nBins = data.size() - 1;
00122   G4double value = 0.;
00123   G4double log_x = std::log10(x);
00124   if (x < points[0])
00125     {
00126       value = 0.;
00127     }
00128   else if (bin < nBins)
00129     {
00130       G4double log_e1 = log_points[bin];
00131       G4double log_e2 = log_points[bin+1];
00132       G4double log_d1 = log_data[bin];
00133       G4double log_d2 = log_data[bin+1];
00134       
00135       //G4cout << "x = " << x << " , logx = " << log_x  << " , bin = " << bin << G4endl; 
00136       //G4cout << "e1 = " << points[bin] << " d1 = " << data[bin] << G4endl;
00137       //G4cout << "e2 = " << points[bin+1] << " d2 = " << data[bin+1] << G4endl;
00138       //G4cout << "loge1 = " << log_e1 << " logd1 = " << log_d1 << G4endl;
00139       //G4cout << "loge2 = " << log_e2 << " logd2 = " << log_d2 << G4endl;
00140 
00141 // Values e1, e2, d1 and d2 are the log values of the corresponding
00142 // original energy and data values. Simple linear interpolation performed
00143 // on loagarithmic data should be equivalent to log-log interpolation
00144       value = log_d1 + (log_d2 - log_d1)*(log_x - log_e1)/(log_e2 - log_e1);
00145 
00146 // Delogarithmize to obtain interpolated value
00147       value = std::pow(10.,value);
00148    }
00149  else
00150    {
00151      value = data[nBins];
00152    }
00153   return value;
00154 }

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