include/Geant4/G4CascadeInterpolator.icc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 #ifndef G4CASCADE_INTERPOLATOR_ICC
0027 #define G4CASCADE_INTERPOLATOR_ICC
0028 //
0029 // Author:  Michael Kelsey <kelsey@slac.stanford.edu>
0030 //
0031 // Simple linear interpolation class, more lightweight than
0032 // G4PhysicsVector.  Templated on number of X-axis (usually energy)
0033 // bins, constructor takes a C-array of bin edges as input, and an
0034 // optional flag whether to truncate or extrapolate (the default) values
0035 // beyond the bin boundaries.
0036 //
0037 // The interpolation action returns a simple double: the integer part
0038 // is the bin index, and the fractional part is, obviously, the
0039 // fractional part.
0040 //
0041 // 20100517  M. Kelsey -- Bug fix in interpolate:  If bin position is _exactly_
0042 //      at upper edge (== last + 0.0), just return bin value.
0043 // 20100520  M. Kelsey -- Second bug fix:  Loop in bin search should start at
0044 //      i=1, not i=0 (since i-1 is the key).
0045 // 20100803  M. Kelsey -- Add printBins() function for debugging
0046 // 20101019  M. Kelsey -- CoVerity reports: recursive #include, index overrun
0047 // 20110728  M. Kelsey -- Fix Coverity #20238, recursive #include.
0048 // 20110923  M. Kelsey -- Add optional ostream& argument to printBins()
0049
0050 #include <iomanip>
0051
0052
0053 // Find bin position (index and fraction) using input argument and array
0054
0055 template <int NBINS>
0056 G4double G4CascadeInterpolator<NBINS>::getBin(const G4double x) const {
0057   if (x == lastX) return lastVal;   // Avoid unnecessary work
0058
0059   G4double xindex, xdiff, xbin;
0060
0061   lastX = x;
0062   if (x < xBins[0]) {           // Handle boundaries first
0063     xindex = 0.;
0064     xbin = xBins[1]-xBins[0];
0065     xdiff = doExtrapolation ? x-xBins[0] : 0.;      // Less than zero
0066   } else if (x >= xBins[last]) {
0067     xindex = last;
0068     xbin = xBins[last]-xBins[last-1];
0069     xdiff = doExtrapolation ? x-xBins[last] : 0.;
0070   } else {              // Assume nBins small; linear search
0071     int i;
0072     for (i=1; i<last && x>xBins[i]; i++) {;}    // Stops when x within bin i-1
0073     xindex = i-1;
0074     xbin = xBins[i] - xBins[i-1];
0075     xdiff = x - xBins[i-1];
0076   }
0077
0078 #ifdef G4CASCADE_DEBUG_SAMPLER
0079   G4cout << " G4CascadeInterpolator<" << NBINS << ">: x=" << x
0080      << " index=" << xindex << " fraction=" << xdiff/xbin << G4endl;
0081 #endif
0082
0083   return (lastVal = xindex + xdiff/xbin);   // Save return value for later
0084 }
0085
0086
0087 // Apply interpolation of input argument to user array
0088
0089 template <int NBINS>
0090 G4double G4CascadeInterpolator<NBINS>::
0091 interpolate(const G4double x, const G4double (&yb)[nBins]) const {
0092   getBin(x);
0093   return interpolate(yb);
0094 }
0095
0096 // Apply last found interpolation to user array
0097
0098 template <int NBINS>
0099 G4double G4CascadeInterpolator<NBINS>::
0100 interpolate(const G4double (&yb)[nBins]) const {
0101   // Treat boundary extrapolations specially, otherwise just truncate
0102   G4int i = (lastVal<0) ? 0 : (lastVal>G4double(last)) ? last-1 : G4int(lastVal);
0103   G4double frac = lastVal - G4double(i);    // May be <0 or >1 if extrapolating
0104
0105   // Special case:  if exactly on upper bin edge, just return value
0106   return (i==last) ? yb[last] : (yb[i] + frac*(yb[i+1]-yb[i]));
0107 }
0108
0109
0110 // Print bin edges for debugging
0111
0112 template <int NBINS>
0113 void G4CascadeInterpolator<NBINS>::printBins(std::ostream& os) const {
0114   os << " G4CascadeInterpolator<" << NBINS << "> : " << G4endl;
0115   for (G4int k=0; k<NBINS; k++) {
0116     os << " " << std::setw(6) << xBins[k];
0117     if ((k+1)%10 == 0) os << G4endl;
0118   }
0119   os << G4endl;
0120 }
0121
0122 #endif  /* G4CASCADE_INTERPOLATOR_ICC */