00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #define INCLXX_IN_GEANT4_MODE 1
00034
00035 #include "globals.hh"
00036
00044
00045 #include <algorithm>
00046 #include <functional>
00047 #include "G4INCLInverseInterpolationTable.hh"
00048
00049 namespace G4INCL {
00050
00051 InverseInterpolationTable::InverseInterpolationTable(IFunction1D const &f, const unsigned int nNodes) {
00052
00053
00054 const G4double x0 = f.getXMinimum();
00055 const G4double x1 = f.getXMaximum();
00056
00057
00058 G4double last = f(x0);
00059 InterpolationNode firstNode(last, x0, 0.);
00060 nodes.push_back(firstNode);
00061 G4int skippedNodes = 0;
00062 for(unsigned i = 1; i < nNodes; i++) {
00063 const G4double xi = x0 + i*(x1-x0)/((G4double)(nNodes-1));
00064
00065
00066 const G4double value = f(xi);
00067 if(value <= last) {
00068 ++skippedNodes;
00069 continue;
00070 }
00071 InterpolationNode node(value, xi, 0.);
00072 nodes.push_back(node);
00073 last = value;
00074 }
00075
00076
00077
00078
00079 initDerivatives();
00080 setFunctionDomain();
00081 }
00082
00083 InverseInterpolationTable::InverseInterpolationTable(std::vector<G4double> const &x, std::vector<G4double> const &y) {
00084
00085
00086
00087
00088
00089 for(unsigned i = 0; i < x.size(); ++i)
00090 nodes.push_back(InterpolationNode(x.at(i), y.at(i), 0.));
00091
00092 initDerivatives();
00093 setFunctionDomain();
00094 }
00095
00096 void InverseInterpolationTable::initDerivatives() {
00097 for(unsigned i = 0; i < nodes.size()-1; i++) {
00098 if((nodes.at(i+1).getX() - nodes.at(i).getX()) == 0.0)
00099 nodes[i].setYPrime(0.0);
00100 else
00101 nodes[i].setYPrime((nodes.at(i+1).getY() - nodes.at(i).getY())/(nodes.at(i+1).getX() - nodes.at(i).getX()));
00102 }
00103 nodes.back().setYPrime(nodes.at(nodes.size()-2).getYPrime());
00104 }
00105
00106 void InverseInterpolationTable::setFunctionDomain() {
00107
00108 if(nodes.front()>nodes.back()) {
00109 xMin = nodes.back().getX();
00110 xMax = nodes.front().getX();
00111 } else {
00112 xMin = nodes.front().getX();
00113 xMax = nodes.back().getX();
00114 }
00115 }
00116
00117 G4double InverseInterpolationTable::operator()(const G4double x) const {
00118
00119 std::vector<InterpolationNode>::const_iterator iter =
00120 std::lower_bound(nodes.begin(), nodes.end(), x);
00121
00122 if(iter==nodes.begin())
00123 return nodes.front().getY();
00124
00125 if(iter==nodes.end())
00126 return nodes.back().getY();
00127
00128 std::vector<InterpolationNode>::const_iterator previousIter = iter - 1;
00129 const G4double dx = x - previousIter->getX();
00130 return previousIter->getY() + previousIter->getYPrime()*dx;
00131 }
00132
00133 std::string InverseInterpolationTable::print() const {
00134 std::string message;
00135 for(std::vector<InterpolationNode>::const_iterator n=nodes.begin(); n!=nodes.end(); ++n)
00136 message += n->print();
00137 return message;
00138 }
00139
00140 }