Geant4  10.01.p03
G4CascadeInterpolator.icc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
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. *
10 // * *
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. *
17 // * *
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 // ********************************************************************
25 //
26 #ifndef G4CASCADE_INTERPOLATOR_ICC
27 #define G4CASCADE_INTERPOLATOR_ICC
28 // $Id: G4CascadeInterpolator.icc 71890 2013-06-27 22:14:13Z mkelsey $
29 //
30 // Author: Michael Kelsey <kelsey@slac.stanford.edu>
31 //
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.
37 //
38 // The interpolation action returns a simple double: the integer part
39 // is the bin index, and the fractional part is, obviously, the
40 // fractional part.
41 //
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()
50 
51 #include <iomanip>
52 
53 
54 // Find bin position (index and fraction) using input argument and array
55 
56 template <int NBINS>
57 G4double G4CascadeInterpolator<NBINS>::getBin(const G4double x) const {
58  if (x == lastX) return lastVal; // Avoid unnecessary work
59 
60  G4double xindex, xdiff, xbin;
61 
62  lastX = x;
63  if (x < xBins[0]) { // Handle boundaries first
64  xindex = 0.;
65  xbin = xBins[1]-xBins[0];
66  xdiff = doExtrapolation ? x-xBins[0] : 0.; // Less than zero
67  } else if (x >= xBins[last]) {
68  xindex = last;
69  xbin = xBins[last]-xBins[last-1];
70  xdiff = doExtrapolation ? x-xBins[last] : 0.;
71  } else { // Assume nBins small; linear search
72  int i;
73  for (i=1; i<last && x>xBins[i]; i++) {;} // Stops when x within bin i-1
74  xindex = i-1;
75  xbin = xBins[i] - xBins[i-1];
76  xdiff = x - xBins[i-1];
77  }
78 
79 #ifdef G4CASCADE_DEBUG_SAMPLER
80  G4cout << " G4CascadeInterpolator<" << NBINS << ">: x=" << x
81  << " index=" << xindex << " fraction=" << xdiff/xbin << G4endl;
82 #endif
83 
84  return (lastVal = xindex + xdiff/xbin); // Save return value for later
85 }
86 
87 
88 // Apply interpolation of input argument to user array
89 
90 template <int NBINS>
91 G4double G4CascadeInterpolator<NBINS>::
92 interpolate(const G4double x, const G4double (&yb)[nBins]) const {
93  getBin(x);
94  return interpolate(yb);
95 }
96 
97 // Apply last found interpolation to user array
98 
99 template <int NBINS>
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
105 
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]));
108 }
109 
110 
111 // Print bin edges for debugging
112 
113 template <int NBINS>
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;
119  }
120  os << G4endl;
121 }
122 
123 #endif /* G4CASCADE_INTERPOLATOR_ICC */