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 #include <algorithm>
00045 #include <cmath>
00046 #include <cstdlib>
00047 #include "G4INCLIFunction1D.hh"
00048 #include "G4INCLLogger.hh"
00049 #include "G4INCLInverseInterpolationTable.hh"
00050
00051 namespace G4INCL {
00052
00053 const G4double IFunction1D::integrationCoefficients[] = {
00054 2.*95.0/288.0,
00055 317.0/240.0,
00056 23.0/30.0,
00057 793.0/720.0,
00058 157.0/160.0,
00059 157.0/160.0,
00060 793.0/720.0,
00061 23.0/30.0,
00062 317.0/240.0,
00063 };
00064
00065 G4double IFunction1D::integrate(const G4double x0, const G4double x1, const G4double step) const {
00066 G4double xi = std::max(x0, xMin);
00067 G4double xa = std::min(x1, xMax);
00068 G4double sign;
00069
00070 if(x1 <= x0) {
00071 sign = -1.0;
00072 std::swap(xi, xa);
00073 } else
00074 sign = 1.0;
00075
00076 const G4double interval = xa - xi;
00077
00078 G4int nIntervals;
00079 if(step<0.) {
00080 nIntervals = 45;
00081 } else {
00082 nIntervals = G4int(interval/step);
00083
00084
00085 G4int remainder = nIntervals % 9;
00086 if (remainder != 0)
00087 nIntervals += 9 - remainder;
00088
00089 nIntervals = std::max(nIntervals, 9);
00090 }
00091
00092 const G4double dx = interval/nIntervals;
00093 G4double result = (operator()(xi) + operator()(xa)) * integrationCoefficients[0]/2;
00094 for(G4int j = 1; j<nIntervals; ++j) {
00095 const G4double x = xi + interval*G4double(j)/G4double(nIntervals);
00096 const unsigned index = j%9;
00097 result += operator()(x) * integrationCoefficients[index];
00098 }
00099
00100 return result*dx*sign;
00101
00102 }
00103
00104 IFunction1D *IFunction1D::primitive() const {
00105 class Primitive : public IFunction1D {
00106 public:
00107 Primitive(IFunction1D const * const f) :
00108 IFunction1D(f->getXMinimum(), f->getXMaximum()),
00109 theFunction(f)
00110 {}
00111
00112 G4double operator()(const G4double x) const {
00113 return theFunction->integrate(xMin,x);
00114 }
00115 private:
00116 IFunction1D const * const theFunction;
00117 } *thePrimitive = new Primitive(this);
00118
00119 return thePrimitive;
00120 }
00121
00122 InverseInterpolationTable *IFunction1D::inverseCDFTable(const G4int nNodes) const {
00123 class InverseCDF : public IFunction1D {
00124 public:
00125 InverseCDF(IFunction1D const * const f) :
00126 IFunction1D(f->getXMinimum(), f->getXMaximum()),
00127 theFunction(f),
00128 normalisation(1./theFunction->integrate(xMin,xMax))
00129 {}
00130
00131 G4double operator()(const G4double x) const {
00132 return std::min(1., normalisation * theFunction->integrate(xMin,x));
00133 }
00134 private:
00135 IFunction1D const * const theFunction;
00136 const G4double normalisation;
00137 } *theInverseCDF = new InverseCDF(this);
00138
00139 InverseInterpolationTable *theTable = new InverseInterpolationTable(*theInverseCDF, nNodes);
00140 delete theInverseCDF;
00141 return theTable;
00142 }
00143
00144 }
00145