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 #ifndef G4CASCADE_INTERPOLATOR_ICC 00027 #define G4CASCADE_INTERPOLATOR_ICC 00028 // $Id$ 00029 // 00030 // Author: Michael Kelsey <kelsey@slac.stanford.edu> 00031 // 00032 // Simple linear interpolation class, more lightweight than 00033 // G4PhysicsVector. Templated on number of X-axis (usually energy) 00034 // bins, constructor takes a C-array of bin edges as input, and an 00035 // optional flag whether to truncate or extrapolate (the default) values 00036 // beyond the bin boundaries. 00037 // 00038 // The interpolation action returns a simple double: the integer part 00039 // is the bin index, and the fractional part is, obviously, the 00040 // fractional part. 00041 // 00042 // 20100517 M. Kelsey -- Bug fix in interpolate: If bin position is _exactly_ 00043 // at upper edge (== last + 0.0), just return bin value. 00044 // 20100520 M. Kelsey -- Second bug fix: Loop in bin search should start at 00045 // i=1, not i=0 (since i-1 is the key). 00046 // 20100803 M. Kelsey -- Add printBins() function for debugging 00047 // 20101019 M. Kelsey -- CoVerity reports: recursive #include, index overrun 00048 // 20110728 M. Kelsey -- Fix Coverity #20238, recursive #include. 00049 // 20110923 M. Kelsey -- Add optional ostream& argument to printBins() 00050 00051 #include <iomanip> 00052 00053 00054 // Find bin position (index and fraction) using input argument and array 00055 00056 template <int NBINS> 00057 G4double G4CascadeInterpolator<NBINS>::getBin(const G4double x) const { 00058 if (x == lastX) return lastVal; // Avoid unnecessary work 00059 00060 G4double xindex, xdiff, xbin; 00061 00062 lastX = x; 00063 if (x < xBins[0]) { // Handle boundaries first 00064 xindex = 0.; 00065 xbin = xBins[1]-xBins[0]; 00066 xdiff = doExtrapolation ? x-xBins[0] : 0.; // Less than zero 00067 } else if (x >= xBins[last]) { 00068 xindex = last; 00069 xbin = xBins[last]-xBins[last-1]; 00070 xdiff = doExtrapolation ? x-xBins[last] : 0.; 00071 } else { // Assume nBins small; linear search 00072 int i; 00073 for (i=1; i<last && x>xBins[i]; i++) {;} // Stops when x within bin i-1 00074 xindex = i-1; 00075 xbin = xBins[i] - xBins[i-1]; 00076 xdiff = x - xBins[i-1]; 00077 } 00078 00079 #ifdef G4CASCADE_DEBUG_SAMPLER 00080 G4cout << " G4CascadeInterpolator<" << NBINS << ">: x=" << x 00081 << " index=" << xindex << " fraction=" << xdiff/xbin << G4endl; 00082 #endif 00083 00084 return (lastVal = xindex + xdiff/xbin); // Save return value for later 00085 } 00086 00087 00088 // Apply interpolation of input argument to user array 00089 00090 template <int NBINS> 00091 G4double G4CascadeInterpolator<NBINS>:: 00092 interpolate(const G4double x, const G4double (&yb)[nBins]) const { 00093 getBin(x); 00094 return interpolate(yb); 00095 } 00096 00097 // Apply last found interpolation to user array 00098 00099 template <int NBINS> 00100 G4double G4CascadeInterpolator<NBINS>:: 00101 interpolate(const G4double (&yb)[nBins]) const { 00102 // Treat boundary extrapolations specially, otherwise just truncate 00103 G4int i = (lastVal<0) ? 0 : (lastVal>last) ? last-1 : G4int(lastVal); 00104 G4double frac = lastVal - G4double(i); // May be <0 or >1 if extrapolating 00105 00106 // Special case: if exactly on upper bin edge, just return value 00107 return (i==last) ? yb[last] : (yb[i] + frac*(yb[i+1]-yb[i])); 00108 } 00109 00110 00111 // Print bin edges for debugging 00112 00113 template <int NBINS> 00114 void G4CascadeInterpolator<NBINS>::printBins(std::ostream& os) const { 00115 os << " G4CascadeInterpolator<" << NBINS << "> : " << G4endl; 00116 for (G4int k=0; k<NBINS; k++) { 00117 os << " " << std::setw(5) << xBins[k]; 00118 if ((k+1)%12 == 0) os << G4endl; 00119 } 00120 os << G4endl; 00121 } 00122 00123 #endif /* G4CASCADE_INTERPOLATOR_ICC */