2 // ********************************************************************
3 // * License and Disclaimer *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
26 #ifndef G4CASCADE_INTERPOLATOR_ICC
27 #define G4CASCADE_INTERPOLATOR_ICC
28 // $Id: G4CascadeInterpolator.icc 71890 2013-06-27 22:14:13Z mkelsey $
30 // Author: Michael Kelsey <kelsey@slac.stanford.edu>
32 // Simple linear interpolation class, more lightweight than
33 // G4PhysicsVector. Templated on number of X-axis (usually energy)
34 // bins, constructor takes a C-array of bin edges as input, and an
35 // optional flag whether to truncate or extrapolate (the default) values
36 // beyond the bin boundaries.
38 // The interpolation action returns a simple double: the integer part
39 // is the bin index, and the fractional part is, obviously, the
42 // 20100517 M. Kelsey -- Bug fix in interpolate: If bin position is _exactly_
43 // at upper edge (== last + 0.0), just return bin value.
44 // 20100520 M. Kelsey -- Second bug fix: Loop in bin search should start at
45 // i=1, not i=0 (since i-1 is the key).
46 // 20100803 M. Kelsey -- Add printBins() function for debugging
47 // 20101019 M. Kelsey -- CoVerity reports: recursive #include, index overrun
48 // 20110728 M. Kelsey -- Fix Coverity #20238, recursive #include.
49 // 20110923 M. Kelsey -- Add optional ostream& argument to printBins()
54 // Find bin position (index and fraction) using input argument and array
57 G4double G4CascadeInterpolator<NBINS>::getBin(const G4double x) const {
58 if (x == lastX) return lastVal; // Avoid unnecessary work
60 G4double xindex, xdiff, xbin;
63 if (x < xBins[0]) { // Handle boundaries first
65 xbin = xBins[1]-xBins[0];
66 xdiff = doExtrapolation ? x-xBins[0] : 0.; // Less than zero
67 } else if (x >= xBins[last]) {
69 xbin = xBins[last]-xBins[last-1];
70 xdiff = doExtrapolation ? x-xBins[last] : 0.;
71 } else { // Assume nBins small; linear search
73 for (i=1; i<last && x>xBins[i]; i++) {;} // Stops when x within bin i-1
75 xbin = xBins[i] - xBins[i-1];
76 xdiff = x - xBins[i-1];
79 #ifdef G4CASCADE_DEBUG_SAMPLER
80 G4cout << " G4CascadeInterpolator<" << NBINS << ">: x=" << x
81 << " index=" << xindex << " fraction=" << xdiff/xbin << G4endl;
84 return (lastVal = xindex + xdiff/xbin); // Save return value for later
88 // Apply interpolation of input argument to user array
91 G4double G4CascadeInterpolator<NBINS>::
92 interpolate(const G4double x, const G4double (&yb)[nBins]) const {
94 return interpolate(yb);
97 // Apply last found interpolation to user array
100 G4double G4CascadeInterpolator<NBINS>::
101 interpolate(const G4double (&yb)[nBins]) const {
102 // Treat boundary extrapolations specially, otherwise just truncate
103 G4int i = (lastVal<0) ? 0 : (lastVal>last) ? last-1 : G4int(lastVal);
104 G4double frac = lastVal - G4double(i); // May be <0 or >1 if extrapolating
106 // Special case: if exactly on upper bin edge, just return value
107 return (i==last) ? yb[last] : (yb[i] + frac*(yb[i+1]-yb[i]));
111 // Print bin edges for debugging
114 void G4CascadeInterpolator<NBINS>::printBins(std::ostream& os) const {
115 os << " G4CascadeInterpolator<" << NBINS << "> : " << G4endl;
116 for (G4int k=0; k<NBINS; k++) {
117 os << " " << std::setw(6) << xBins[k];
118 if ((k+1)%10 == 0) os << G4endl;
123 #endif /* G4CASCADE_INTERPOLATOR_ICC */